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ABSTRACT 

We  discuss  the  three  fundamental  issues  of  a  computational  approach  in 
structure  prediction  by  energy  minimization  and  examine  them  for  the  nucleic 
acid  component  deoxyribose.  Predicting  the  conformation  of  deoxyribose  is 
important  not  only  because  of  the  molecule's  central  conformational  role  in 
the  nucleotide  backbone,  but  also  because  energetic  and  geometric 
discrepancies  from  experimental  data  have  exposed  some  underlying 
uncertainties  in  potential  energy  calculations.  The  three  fundamental  issues 
examined  here  are:  i)  choice  of  coordinate  system  to  represent  the  molecular 
conformation,  ii)  construction  of  the  potential  energy  function,  and  iii) 
choice  of  the  minimization  technique.  Our  approach  consists  of  the  following 
combination.  First,  the  molecular  conformation  is  represented  in  cartesian 
coordinate  space  with  the  full  set  of  degrees  of  freedom.  This  provides  an 
opportunity  for  comparison  with  the  pseudorotation  approximation.  Second, 
the  potential  energy  function  is  constructed  so  that  all  the  interactions  other 
than  the  non-bonded  terms  are  represented  by  polynomials  of  the  coordinate 
variables.  Third,  two  powerful  Newton  methods  that  are  globally  and 
quadratically  convergent  are  implemented:  Gill  and  Murray's  .Modified 
Newton  method  and  a  Truncated  Newton  method,  specifically  developed  for 
potential  energy  minimization.  These  strategies  have  produced  the  two 
experimentally-observed  structures  of  deoxyribose  with  geometric  data  (bond 
angles  and  dihedral  angles)  in  very  good  agreement  with  experiment.  .More 
generally,  the  application  of  these  modeling  and  minimization  techniques  to 


potential  energy  investigations  is  promising.  The  use  ol  cartesian  variables 
and  polynomial  representation  of  bond  length,  bond  angle  and  torsional 
potentials  promotes  efficient  second-derivative  computation  and,  hence, 
application  of  Newton  methods.  The  truncated  Newton,  in  particular,  is 
ideally  suited  for  potential  energy  minimization  not  only  because  the  storage 
and  computational  requirements  of  Newton  methods  are  made  manageable, 
but. also  because  it  contains  an  important  algorithmic  adaptive  feature:  the 
minimization  search  is  diverted  from  regions  where  the  function  is  non- 
conve.x  and  is  directed  quickly  toward  physically  interesting  regions. 
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1.     INTRODUCTION 

Theoretical  approaches  for  predicting  three-dimensional  structures  of  nucleic  acids  and 
proteins  are  now  recognized  as  powerful  tools  for  revealing  details  of  molecular  conformation, 
motion,  and  associated  biological  functions.  As  a  recent  article  in  Science  states.  "It  is  clear  that 
theoretical  chemistry  has  entered  a  new  stage  ...  with  the  goal  of  being  no  less  than  full  partner 
with  experiment"  [1]. 

Indeed,  theoretical  calculations  using  semi-empirical  potential  energy  functions  can 
complement  information  obtained  from  experimental  techniques.  >uch  as  X-Ray  crystallography 
and  spectroscopic  methods  regarding  structural  features  of  a  molecule,  and  can  make  reliable 
predictions  [1--*].  The  fundamental  idea  behind  potential  energy  calculations  is  to  con^truct  a 
function  that  represents  the  free  energy  associated  with  a  specific  molecular  conformation;  the 
conformation  corresponding  to  the  minimum  free  energy  can  thenjtte.  appiroximated  by  potential 
energy  minimization.  The  challenge  and  power  of  this  approach  lies-in»  its-ability  to  incorporate  all 
experimental  data  known  from  crystal  structures  of  small  molecules  with  the  goal  o(  predicting 
structures  of  large  biological  systems.  Understanding  the  correlation  between  conformation  and 
biological  function  is  our  ultimate  goal. 

Two    basic    ingredients  of    potential    energy    studies    make    these    investigations    difficult: 

construction    of    the    energy  function    and    energy    minimization.    With    increasing    accuracy    of 

experimental  data  and  better  computing  resources,  these  difficulties,  however,  may  be  overcome. 

Nonetheless,     understanding  the    basic    computational     problems     is     important     for     improved 
theoretical  treatments. 

[n  this  paper,  we  discuss  the  three  fundamental  issues  of  a  computational  approach  —  choice 
of  coordinate  system  to  represent  molecular  conformation,  constructum  ot  the  potential  energy 
function,  and  choice  of  minimization  technique.  We  then  present  the  particular  combination  that 
we  have  devised  and  applied  to  a  deoxyribose  model.  The  preferred  sugar  conformations  of 
C2'  — endo  and  C3"  — endo  (see  Fisures   1.2)  are  venerated  with  aeometries  as  observed  for  manv 


sugar  t'ragmLMits  of  crystallographicall\-determined  nucleoNides  and  nucleotides  [5-9j. 

Deo\>rib<')se  is  particularl_\'  suitable  as  an  application  model  since  it  is  small  enough  _\et 
complicated  and  biologically  important.  Linked  to  both  the  phosphate  and  base  group  of  a 
nucleotide,  the  furanose  ring  conformation  strongly  influences  the  overall  conformation  of  a 
nucleic  acid.  Mathematically,  the  problem  is  challenging  because  the  ring  atom  coordinates  are 
difficult  to  generate  due  to  constraints  imposed  by  ring  closure  [10-161.  Manv  earl>  attempts  to 
reproduce  the  observed  puckering  preferences  of  nucleic  acid  sugars  wiih  associated  endocyclic 
bond  angle  values  by  potential  energ\  methods  uere  imperfect  [9].  These  discrepencies  were 
attributed  to  the  transfer  of  energy  parameters  from  general  chemical  sequences  in  small 
molecules  to  analogous  chemical  groups  in  the  furanose  ring^  Olson  has  reconciled  these 
differences  by  incorporating  directly  experimental  values  of  the  bond  angles  that  acct^mpany  the 
puckering  and  by  improving  the  energy  function  with  a  'gauche  potential  [I7j.  While  more  recent 
minimization  or  dynamics  studies  [18-20]  have  been  more  successful  at  reproducing  qtialirative 
aspects  of  pseudorotation  (regions  of  energy  minima  and  maxima,  values  of  energy  barriers  i. 
there  are  still  some  discrepancies  from  experiment  in  the  exact  location  of  the  minima  and 
maxima,  and.  when  reported,  in  the  bond  angle  deviations  that  occur  with  puckering.  .Accurate 
generation  of  structures  and  energies  for  the  furanose  ring  by  full  cartesian  energy  minimization  is 
important  in  order  to  compute  nucleic  acid  structures  by  minimization  in  conformational  space 
with  all  degrees  of  freedom  rather  than  in  dihedral  angle  space  [see.  for  example.  21-23]. 

Our  computational  approach  consists  of  the  following  three  components.  First,  the  potential 
energy  is  represented  in  cartesian  coordinate  space  with  the  full  set  of  degrees  o[  freedom 
associated  u  ith  the  molecular  conformation.  Thus,  all  bond  lengths  and  bond  angles  as  well  as 
dihedral  angles  are  allowed  to  vary.  When  properly  parameten/.ed.  this  approach  provides  a 
realistic  view  of  the  relaxed  molecular  conformation.  .Second,  the  potential  energy  function  is 
cimstructed  so  that  all  the  interactions  other  than  the  non-bonded  terms  are  represented  by 
polynnmials  of  the  coordinate  variables.  This  formulation  allov^s  direct  differentiation  vMth  respect 
to  the  cartesian   variables  and   facilitates   manipulation   and   differentiation   of  the  energy  terms. 


Third,  two  powerful  Newton  methods  that  are  globally  and  quadratically  convergent  [2-1  are 
unpleinented:  Gill  and  Murray's  Modified  Newton  method  [2>.26j  and  a  Truncated  Neuton 
method  specifically  developed  for  potential  energy  minimization.  Full  details  of  the  minimization 
algorithms  are  provided  in  the  accompanying  paper  [68]. 

In  the  ne.xt  section,  we  describe  the  fundamental  issues  of  computational  approaches  while 
explaining  the  directions  we  have  chosen.  Section  3  discusses  the  methods:  geometric  construction, 
energy  formulation  and  parameterization,  and  organization  for  minimization.  In  Section  4.  ue 
examine  in  detail  the  influence  of  different  potential  energy  sets  on  the  structures  obtained  and  on 
the  C3'  — endo'C2"-endo  energy  difference.  In  particular,  ue  can  study  the  contribution  of  each 
energy  potential  from  Energy  vs.  P  (the  pseudorotation  parameter)  profiles.  These  curves  also 
provide  an  opportunity  to  investigate  the  approximation  made  in  the  pseudorotation  description. 
In  Section  5  we  summarize  our  overall  conclusions.  S- 
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2.     COMPUTATIONAL  APPROACHES 

In  anv  computational  approach  to  the  determination  of  molecular  >tructure  by  energy 
minimization,  three  basic  decisions  must  be  made: 

1)  The  degrees  of  freedom  used  to  describe  the  molecular  structure. 

2)  The  analytic  form  of  the  energy  function,  and 

3)  The  choice  of  a  minimization  algorithm. 

The  structural  outcome  and  hence  the  biological  implication  of  any  calculation  is  highly  dependent 
on  the  combination  of  choices  taken.  Given  the  ^ame  starting  point,  different  minimization 
algorithms  may  predict  different  local  minima,  and  different  structures  may  result  froin  different 
potential  energy  models.    We  will  describe  these  three  issues  in  turn. 

2.1     Degrees  of  Freedom 

The  conformation  of  a  molecule  is  described  by  a  list  of  numbers  that  specifies  the  relative 
positions  of  the  atoms  m  space.  By  definition,  the  conformation  is  unchanged  uhen  the 
molecule  as  a  whole  is  subjected  to  rigid-body  motion  (translation  or  rotation).  If  the  molecule 
contains  .V  atoms.  3iV-6  numbers  are  required  to  specify  the  conformation.  These  numbers 
may  be  chosen  in  different  uays.  For  example,  one  might  use  the  3.V  cartesian  coordinates  o\  the 
atoms  and  force  6  of  these  numbers  to  be  zero  by  putting  a  particular  atom  at  the  origin,  a  second 
atom  along  the  v  a.xis.  and  a  third  atom  in  the  v  .  \  plane.  .Alternatively .  the  contormation  may  be 
.specified  as  some  combinatii)n  of  bond  lengths,  bond  angles,  and  dihedral  angles  |27]. 

The  internal  energy  of  a  molecule  is  some  function  of  its  conformation  i  the  energy  is 
unchanged  by  translation  or  by  rotation  oi  the  molecule  as  a  v\hiile).  Thus  the  task  of  energy 
minimization  must  be  carried  out  in  conformation  space,  which  has  3A'-6  dimensions  (degrees  of 


freedom).  Since  this  number  of  dimensions  can  be  very  large  for  biological  macromolecules.  it  is 
temptmg  to  reduce  the  number  of  degrees  of  freedom  by  makmg  use  ot  certain  a  priori 
information  concerning  the  molecule  in  question.  One  might,  for  example.  as>ume  that  the  bond 
lengths  and  bond  angles  were  rigid  and  known  in  advance.  Then  the  conformation  of  the  molecule 
(and  hence  the  energy)  would  be  a  function  of  the  dihedral  angles  alone.  .Another  possibility 
would  be  to  assume  that  the  molecule  could  be  deformed  only  along  >ome  particular  path  in 
conformation  space.  The  use  of  the  pseudorotation  model  is  an  example  of  the  latter  approach. 

The  concept  of  pseudorotation  [5. 1  1 .  12. 14..28.29]  restricts  the  energetic  pathway  that  the 
five-membered  sugar  ring  follows  to  a  wave-like  motion  from  a  chosen  mean  plane  detined  by  the 
five  ring  atoms.  The  conformation  of  these  skeletal  atoms  is  described  by  only  two  instead  k)\  the 
full  9  degrees  of  freedom  (3.V-6  ,  ,V  =  5)  a^soclated  u  ith  the  molecular  conformation.  These 
coordinates  can  be  calculated  in  the  Cremer  and  Pople  formalism  [U]  by  expressing  the 
r  — coordinates  of  the  atoms  as  periodic  displacements  from  a  chosen  mean  plane  in  terms  of  phase 
amplitude  and  phase  shift  (t/.'^}: 

.-^  =  (2/5)'-c/ cos  (  ^  -r  -^  (;-2)  )  .  ;  =  0.1.2,3.4.  (1) 

.Alternatively,  the  coordinates  can  be  constructed  in  the  .Altona  and  Sundaralingam  description  [5] 
from  the  ring's  dihedral  angles,  expressed  as  periodic  variations  in  terms  of  amplitude  and  shift 
("ma.x-  P)  '^ee  Figure  3): 

T;  =  T^^,  cos(  P  +  ^(;-2)  )  .  7  =  0.1.2.3.4.  (2) 

.Although  application  of  this  concept  to  nucleic  acids  is  conceptually  Mmple  and  can  simplify 
formulation  of  ring  geometry,  it  is  necessarily  approximate.  First,  it  ma\  produce  anomalies  in  the 
overall  nucleic  acid  structure.  Second,  the  pseudorotation  formulation  introduces  mathematical 
difficulties  in  energy  representation  fas  a  function  of  the  2  pseudorotation  parameters)  and  energy 
differentiation  (with  respect  to  the  conformational  variables). 


Thus,  while  it  is  generally  an  advantage  to  work  u  iih  t'ewer  degrees  of  freedom,  there  are 
several  problems  associated  with  the  use  of  constraints: 

1)  The  constraints  may  be  unrealistic.  In  a  real  molecule,  bond  lengths  and  particularly  bond 
ansles  can  fluctuate  in  response  to  molecular  forces  that  vary  v\ith  conformation.  Such  >mall 
changes  may  produce  large  deformations  of  the  molecule  as  a  \\hole  [9]. 

2)  The  constraints  may  be  inconsistent!  This  is  particularly  true  in  the  case  of  five-atom  rings, 
for  which  the  usual  constraints  on  bond  lengths  and  bond  angles  are  inconsistent  with  ring 
closure  [10]. 


3)  When  constraints  are  used  to  reduce  the  number  of  independent  \ariables,  the  energy  may  be 
a  very  complicated  function  of  the  variables  that  remain.  In  particular,  the  non-bonded 
interactions  (Coulomb  and  Van  der  Waalsl  are  most  easily  expressed  as  a  function  of  the 
cartesian  coordinates  of  the  atoms.  If.  for  example,  the  independent  variables  are  taken  as  a 
collection  of  dihedral  angles,  then  the  cartesian  coordinates  must  be  computed  from  the 
dihedral  angles  at  every  stage  of  the  minimization  process.  Worse,  if  the  minimization 
method  uses  derivatives,  the  corresponding  derivatives  of  the  cartesian  cocudinatcs  uith 
respect  to  the  dihedral  angles  must  be  evaluated  repeatedK  .  The  chain-rule  expressions  for 
the  derivatives  (especially  second-derivatives)  are  tedious  to  derive,  easy  to  get  wrong,  and 
expensive  to  evaluate.  Symbolic  computation  (e.g.  Macsyma  [30])  avoids  the  tedium  and  the 
potential  errors,  but  not  the  expense  of  repeatedly  evaluating  the  resulting  expressions. 

For  all  ot  the  above  reasons,  we  avoid  constraints  and  work  uith  the  cartesian  coordinates  of 
the  atoms  as  the  independent  variables.  Our  energy  function  contains  terms  which  keep  the  bond 
lengths  and  bond  angles  close  to  their  observed  values,  but  the  stiffness  o{'  these  "soft  constraints" 
is  not  infinite,  so  we  avoid  the  pitfalls  outlined  above. 
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2.2     The  Potential  Energy  Function 

In  principle,  a  potential  energy  function  could  be  obtained  by  a  quantum-mechanical 
description  of  the  ground-state  energy  of  the  molecule.  Since  >uch  calculations  are  not  yet  feasible 
for  molecules  as  large  as  proteins  and  nucleic  acids,  ue  resort  to  a  molecular  mechanics  treatment. 
We  consider  the  molecule  as  a  system  of  .V  masses  (atoms)  which  is  deformed  by  molecular 
forces,  thereby  producing  an  energy  change.  The  basic  form  of  this  energy  consists  of  a  sum  of 
non-bonded  interactions,  torsional  potentials,  and  strain  energy: 


^.\o.\ 


BOSDED 
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^sos -BOSDED  represents  the  pairvMse  interactions  and  consists  of  two  contributions:  Van  der 
Waals  interactions,  repulsive  at  short  atomic  separations  and  attractive  at  large  distances:  and 
electrostatic  interactions  between  charged  groups  obeying  Coulomb's  law.  The  Lennard-Jones 
6-12  empirical  potential  was  formed  by  combining  the  leading  term  (  -.4  /  r^  )  o(  the  London 
attraction  potential,  knov\.n  from  quantum  mechanical  theory,  with  a  steep  repulsion  term  that  is 
computationally  convenient.  Interactions  are  generally  considered  for  atom  pairs  in  .S'vg.  the  set  oi 
all  (/.;)  pairs  for  which  ;  <  j  and  atoms  /  and  j  are  separated  by  3  bonds  or  more.  This  avoids 
double  consideration  (in  terms  (3a)  and  (3c))  of  bonded  atoms  and  atoms  involved  in  a  bond 
angle  in  some  calculations,  to  avoid  double  consideration  of  atoms  involved  in  a  dihedral  angle  as 
well.  .')yfl  ci^nsists  of  the  atom  pairs  separated  by  4  bonds  or  more.  Thus.  v\hcn  atoms  separated 
by  3  bonds  e.xactly  are  considered  non  — bonded .  the  contribution  to  the  torsional  potential  must  be 
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recoenized:  this  can  be  done  by  adjusting  the  torsional  parameters  or  scaling  the  1  --  interactions. 
The  electrostatic  potential  energy  between  the  partial  atomic  charges  Q,  is  given  b> 
Coulomb's  law  modified  by  a  dielectric  function  Dir,.]  to  account  for  a  weaker  interaction  in  a 
polarizable  medium  than  in  a  vacuum.  For  convenience,  summation  generally  extends  over  atoms 
defined  as  non-bonded  for  the  Lennard-Jones  potential.  This  involves  an  approximation,  hecau.se 
although  we  attempt  to  avoid  the  contribution  of  atom  pairs  separated  by  1  or  2  chemical  bonds  in 
more  than  one  energy  term,  this  convention  assumes  neutral  interactions  between  these  excluded 
atom  pairs.  ■ 

The    torsional    potential    ^torsio^    accounts     for    interaction     between    atoms     invoKed     in 

I 
internal  rotation .   a   rotation   about   a   bond   connecting   two   chemical   groups    in    a   molecule   and 

described  by  a  dihedral  angle  7  (see  Figure  -i-)-  In  ethane,  for  example,  the  C  — C  bond  provides  a 

rotation  a.xis  for  the  2  methyl  groups.    .Although  the  origin  of  the  barrier  to  internal  rotation  has 

not   been   entirely   resolved,   the    principal    interactions   that   give    rise   to   rotational   barriers   are 

currently  thought  to  be  repulsive  interactions,  caused  by  overlapping  of  bond  orbitals  of  the  two 

rotating  groups  [31].    In  ethane,  the  torsional  energy  is  highest  when  the  two  ineth\l  groups  are 

nearest,  as  in  the  eclipsed  state,  and  lowest  when  the  two  groups  are  optimall>  separated,  as  m  the 

staggered      state.        The       empirical       form       of      the       torsional       potential       is      given       by 

E-TOR  ('^^  —  ~(  1  +COS/ITJ,    where    the    integer    n    denotes    the    periodicity    ot    the    rotational 

barrier    and    V^    is    the    associated    barrier    height.     Tv\ofold    and    threefold    potentials    are    most 
commonly  used. 

^STRM.\  represents  the  bond  stretching  and  angle  bending  energy  as  bond  lengths  and  bond 
angles  (  8,  )  deviate  from  equilibrium  values  (denoted  in  the  energy  formulas  with  bars).  .V^ 
denotes  the  set  of  all  ii.j)  pairs  for  which  i<j  and  the  atoms  are  bonded.  .Since  chemical  bonds 
display  a  very  narrow  range  of  fluctuation  in  length  i  generall\  in  the  order  o\  0.1  .-^  i .  a  harmonic 
or  quadratic  approximation  is  considered  sufficient.  Bond  angle  deviations  are  usuall>  small 
(<  3")   unless  electron  lone  pairs  t\)rm  bond-like  orbitals  (e.2.  water.  B  (  H  - () - H  )  ~   105")  or 
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ring  molecules  impose  closure  constramrs  (e.g..  cycloalkanes.  deoxyribose) .  For  these  larger 
deviations,  suitable  potentials  are  still  bemg  investigated,  but  the  most  commonly  used  form  of  the 
bending  potential  is  also  harmonic. 

.Additional  terms  or  variations  to  the  basic  energy  form  in  (3)  may  be  used  to  account  for 
hydrogen  bonding,  solvent  effects  or  helical  parameters  where  appropriate  [  18. 32. 33. for 
example]. 

The  parameterization  process  for  potential  energy  functions  is  a  difficult  task  Several 
important  decisions  must  be  made  regarding  choices  for  the  functional  form  and  numerical  values 
for  the  parameters.  Even  if  one  is  given  a  specific  energy  form  and  a  set  of  structural  and 
energetic  data  to  reproduce,  the  combination  of  parameters  that  can  be  used  are  endless. 
Unrealistic  choices  for  one  group  ot  parameters  can  be  compensated  for  by  adjustment  of  another. 
Ironically,  the  more  specific  the  force-field  becomes,  that  is,  the  more  individualized  the  treatment 
of  geometric  sequences  (bonds,  bond  angles  and  dihedral  angles),  the  greater  becomes  the 
possibility  of  divergence  from  physical  reality. 

[n  theory,  the  energy  terms  ^hould  have  clear  physical  significance  with  parameters  calibrated 
by  empirical  fitting  of  crystal  data  and  rotational  barriers  of  analogous  small  molecules  However, 
an  approximation  is  inherent  in  the  extension  of  data  from  small  to  large  ^ystems.  .Moreover, 
interaction  with  solvent  and  countenons,  reflected  in  the  experimental  data,  must  be  interpreted 
and  incorporated  in  the  energy  model.  In  summary,  much  freedom  and  manipulation  are  possible 
in  constructing  semi-empirical  energy  surfaces.  Only  if  constructed  and  parameterized  correctly, 
will  the  energy  model  generate    reliable  structural  predictions. 

In  the  case  of  nucleic  acid  sugars,  the  importance  of  parameter  choices  in  the  energy  function 
has  already  been  realized.  Different  potential  energy  models  have  produced  revults  that  are 
cjiialitatively  different  -  pseudorotation  is  "free"  [34|.  or  pseudorotation  is  hindered  [18-20].  This 
is  particularly  possible  when  chosing  equilibrium  values  for  endocyclic  bond  angles  that  are 
inappropriate  [17.34-].  Since  the  unusual  puckering  geometry  and  ring  closure  constraints  produce 
significant  deviations  from  tetrahedral  bond  angle  arrangements,  it  is  not  clear  uhat  equilibrium 
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values  should  be  used  in  the  harmonic  bending  terms,  more  appropriate  for  small  fluctuations. 
Indeed.  Olson's  studies  [17]  were  >uccesstul  at  producmg  energy  minnna  at  the  ideal 
pseudorotation  phase  angles  of  18'^  and  162".  partl_\  because  they  mvolved  no  mmmiizaiion; 
experimentally-determined  valence  angles  were  directly  mcorporated  mto  the  coordinate 
generation  oi  the  ring  atoms  to  give  the  energy  as  a  function  of  a  fixed  puckered  form. 

In  addition  to  the  sensitivity  of  the  results  to  the  parameter  choices,  important  issues  m  the 
modeling  procedure  have  emerged  in  deoxyribose  investigations.  For  example,  the  drbitrar\ 
selections  possible  in  modeling  dihedral  angles  have  been  noted  by  Harvey  and  Prabhakaran  [l'^]. 
Consider  an  ethane  molecule;  nine  rotations  are  associated  uith  the  C-C  bond.  Which  dihedral 
angles  should  we  associate  with  a  particular  torsional  term.'  .And  how  should  we  assure  proper 
positioning  of  the  remaining  atoms'.'  Since  the  torsional  potentials  have  the  greatest  contribution  to 
the  total  energy 'of  deoxyribose.  Harvey  and  Prabhakaran  suggest  that  a  ver\  careful  search  for 
appropriate  torsion  parameters  and  modeled  dihedral  angle  sets  are  required  to  obtain  energy 
minima  at  the  angles  corresponding  to  the  ideal  C2"  — endo  and  C3'-endt)  structures.  Selection  of 
the  torsion  parameters  is  particularly  important  for  deoxyribose.  since  a  combination  ot  both 
threefold  and  twofold  (gauche)  torsional  energies  is  required  to  produce  relative  energies  ot  the 
preferred  C2'-endo  and  C3"  — endo  structures  in  accordance  with  the  observed  ratio  of  these 
forms  for  sugar  fragments  [5-9].  We  will  discuss  the  selection  of  modeled  dihedral  angles  and 
associated  constants  in  Section  3.  .Although  choices  in  the  torsion  parameters  and  modeled 
dihedral  angles  are  important  for  the  furanose  ring,  the  discrepancies  noted  in  [19]  regarding 
location  of  the  minima,  are  undoubdetly  attributed  in  part  to  the  united  atom  approach,  Rxplicit 
inclusion  ot  the  hydrogen  atoms  is  necessary  in  order  to  account  correctly  for  the  eclipsing  strain 
of  the  two  substituents  of  C2'  and  C3"  in  the  East  and  the  \'an  der  Waals  repulsions  between  the 
equatorial  base  (in  parallel  orientatum  to  the  furanose  ring  plane)  and  the  exoc\clic  substituents  of 
C5'  in  the  West. 

In  our  study,  all  heavy  atom  rotations  about  one  bond  (rotations  not  involving  hvdrogen 
atoms)  are  considered  systematically,  and  all  force-field  constants  are  used  as  group  parameters 


-  11  - 

tor  generally-defined  families  of  atomic  sequences.  This  procedure  allows  the  molecule  to  attain 
equilibrium  bond  lengths,  bond  angles,  and  dihedral  angles  as  a  result  of  energy  minimization 
rather  than  targeting  strategies.  It  can  also  help  to  prevent  generation  of  unrealistic  geometries. 

2.3     Minimization 

For  potential  energy  minimization  the  choice  of  method  is  important  because  of  two  inherent 
difficulties:  existence  of  many  local  minima,  all  possibly  of  biological  importance,  and  extensive 
computational  requirements.  The  state-of-the-art  today  is  such  that  for  many  small  problems 
(about  30  variables  or  less)  suitable  algorithms  exist  for  finding  all  local  minima.  For  larger 
problems,  however,  unless  good  initial  approximations  are  provided,  there  is  no  guarantee  of 
solving  them  completely  in  a  finite  number  of  trials.  Many  trials  are  generally  required  to  find 
various  local  minima,  and  finding  the  global  minimum  cannot  be  as.sured. 

The  choice  of  a  minimization  algorithm  must  be  made  by  considering  the  following  four 
features  of  the  problem: 

1)  Form  of  the  objective  function.    Is  it  linear  or  nonlinear'.^  smooth  or  discontinuous?  a  sum  of 
squares?  with  or  without  constraints'.'  with  or  uithout  bounds? 

2)  Size.     For    a    large    number    of   variables    storage    considerations    are    important.    Effective 
techniques  for  small  problems  are  usually  unsuitable  for  large  problems. 

3)  .Availability  of  analytic  derivatives. 

4|      Computer  resources.    The   better  the  storage  and  :^peed  capabilities,  the  more   flexible  the 
choice  of  an  algorithm. 

Potential  energy  functions  are  generally  large,  nonlinear,  and  the  problem  can  be  formulated 
as  unconstrained.    Obtainins  analvtic  first  and  >econd-derivatives  inav  be  difficult  but  is  definitelv 


feasible.  Since  storage  and  cost  considerations  are  directly  related  to  the  comple\it\  of  function 
and  derivative  calculations,  the  choice  of  a  miniinization  algorithm  for  potential  energ_\  functions 
should  be  based  on  the  availability  of  derivatives.  We  classify  minimization  algorithms  into  three 
categories: 

(i)    Non-derivative    methods, 
(ii)    First-derivative  (or  gradient)  methods,  and 
(iii)    Second-derivative  methods. 

Non-derivative  methods  (e.g  PowelTs)  are  generally  easy  to  implement,  but  they  tend  to  be 
inefficient  and  e.xcessively  slow.  The  computational  cost,  dominated  b_\  the  number  ot  function 
evaluations,  can  be  large  for  functions  of  many  variables  and  thus  far  outweigh  the  benefit  of 
avoiding  derivative  computation.  Gradient  methods  (e.g.  nonlinear  Conjugate  Gradient)  are  the 
most  commonly  used  in  potential  energy  minimization  [  lS.32.35r38.  for  example]  because  their 
storage  and  computational  requirements  are  generally  manageable  and  their  convergence 
properties  satisfactory,  although  often  unpredictable. 

Newton-based  or  .Nccond-derivative  methods  are  particularly  attractive  tor  nonlinear 
optimization  problems  because  of  their  locallx  rapid  convergence  properties.  With  appropriate 
global  modifications,  convergence  to  a  local  minimum  can  be  guaranteed  e\en  from  distant  initial 
coordinates.  However,  the  use  of  Newton  methods  has  been  restricted  to  small  problems  or  only 
near  a  solution  following  a  gradient  method  [32.35.3<S.  for  example]  as  a  result  of  the  demanding 
computational  and  storage  requirements  associated  with  calculation  and  manipulation  oi  the 
second-derivative  matrix.  Fortunately,  advances  in  computing  resources  and  algorithms  are 
making  Newton  methods  feasible  and  powerful  for  large-^cale  problem^.  Trmuati'd  Newton 
methods,  in  particular,  can  maintain  the  quadratic  convergence  of  Nev\ton  methods  uhile 
considerably  reducing  the  co>tly  manipulation  and  demanding  .^torage  of  the  analytic  Hes>ian 
matrix. 


In   Appendix   A.   we   briefly   highlight  the   general  concepts  of  minimization   >trategies.  and 
refer  the  reader  for  more  details  to  selected  references  [39-44], 


3.     METHODS 


3.1     Geometry 

We  consider  a  model  tor  the  2"  — deoxyribose  ring  (Figure  1)  v.  ith  the  full  set  of  degrees  of 
freedom  associated  with  the  5  rmg  atoms.  Since  we  are  interested  in  the  structural  behavior  o[ 
this  constituent  v\hen  linked  to  a  base  and  phosphate  groups,  we  consider  the  NH,.  OH.  and  CH, 
attachments  as  atom  groups.  Hvdrogens  linked  to  ring  atoms  are  treated  as  mdi\idual  atoms  and 
not  as  part  of  united  CH  or  CH^  groups.  Individual  treatment  of  these  h\drogens  is  necessary  in 
order  to  represent  accurately  the  non-bonded  energy  associated  with  the  molecule.  Without 
individual  representation  of  the  C2'  hydrogens,  for  example,  no  barrier  will  occur  for  the  non- 
bonded  energy  component  in  the  East  region  of  the  pseudorotation  path.  v\here  one  C2'  hydrogen 
IS  eclipsed  with  the  C3"  hydroxyl  group.  Thus,  our  model  consists  of  thirteen  atomic  positions  in 
cartesian  coordinate  space. 

We  express  bond  lengths  and  cosines  of  bond  angles  and  dihedral  angles  as  polynomials  of 
the  cartesian  variables  using  a  "hierarchy"  of  expressions  as  follows: 

Let  X,  =  (   v,  J  .  V,  2  ■  V,  3  )  .  /  =   1 V    (  .V  is  the  number  of  atoms   )  denote  the  p<-)sition 

vector  of  atom  (  and  r,j  =  x,  -  x^  .  denote  the  distance  vector  from  atom  j  to  /.  For  any  vector 
a  we  denote  the  magnitude  by  {|ai{  and  the  associated  equilibrium  magnitude  by  a.  Similarly,  if 
an  equilibrium  value  is  given  for  a  bond  angle  9  .  ue  denote  it  by  9  .  For  convenience  later  in 
writing  the  potential  energy  function,  we  denote  the  vector  magnitude  ||r,J|  also  as  i\,.  in  non-bold 
type. 

We  then  express  the  cosine  of  a  bond  angle  9,,^,  t\)rmed  by  atoms  i-j-k  as  an  inner  product 
with  equilibrium  bond  lengths: 

(x^-x.)     (X,  — x^, ) 
cos  9,^;.  =  z — z ^~.  (-l-) 

'■■7     ''^7 


Let  -,,1:1  be  the  dihedral  angle   defining  the   rotation   about  bond   j-k  (see   Figure  4).     To 
»nTiplity  the  notation  we  define  the  following  quantities: 


b  =  r^,^  =  Xk-\j 

c  =  r,;.  =  x,-x;, 

9^^  —  angle  between  a  and  b 

9;,^  —  angle  between  b  and  c 

n^i,  —  unit  normal  to  plane  spanned  by  \ectors  a  and  b 

ri/,^  —  unit  normal  to  plane  ^panned  by  \ectors  b  and  c. 


Then  by  definition,  x,^;./  is  the  angle  between  the  normals  n^,^  and  n^^ 


cos  -,ji,,  -  n 


or 


^ab       "ftc 


a  X  b 


cos   7 


ijkl 


b  X  c 


|a|  |b|  sin  9^^       |bi  |c|  sin  9 


be 


The  sign  of  7,^^./  is  determined  by  the  sign  of  the  triple  product  (a  x  bi  ■  c  (see  Figure  4).  To 
obtain  an  eighth-degree  polynomial  in  the  coordinate  variables,  we  replace  bond  lengths  ||a||.  |ib|| 
and  i|c||  and  bond  angles  9^^^  and  9;,^  by  their  equilibrium  values  and  then  simplify  this  expression 
to  a  difference  of  inner  products: 


(a  ■  bi  (b  ■  c)  -  fa  ■  CI  (b  ■  bi 


cos 


tjkl 


a    b  ■    c    sin  9  ,/,  sin  9 


(5) 


he 


This  geometric  representation  is  direct  and  Mmple  for  energ\'  e\'aluatiim  and  differentiation 
and  is  motivated  by  the  physical  "hierarchy"  of  flexibility  -  bond  lengths  are  nearly  constant,  bond 
anales  are  somewhat  variable,  and  dihedral  ancles  are  the  most  flexible  in  their  ranees. 
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3.2    Potential  Energy  Function 

The  conformational  energy  of  the  deoxyribose  model  is  constructed  as  a  sum  of  contributions 
from  a  standard  6-12  Lennard-Jones  function  iEi_j\.  a  coulombic  potential  iEcocL^-  bond  length 
and  bond  angle  strain  terms  iEgQxp  and  f^s.A.vG'-  ^""^  torsional  potentials  with  rotational  barrier 
periodicity  of  two  and  three  (EifoLD  ^"'l  Eyfo^Q): 


^   "   ^U    "•"   ^COLL    ^   ^BOSD    "•"   ^BAXG    +   ^ IFOLD    "*"   ^3 


}FOLD 


(6) 


'■J   i  ^sa 


+ 


i: 


I  6a) 


Q,Q, 


e  e       r, 


(  3  =  0.1.     £  =  4.0 


(6b) 


I.J  i  s. 


t6c 


^B,\.\G   -    S 


52,   [  cos  9,  -  cos  9,  ) 


(6d) 


E^ 


FOLD 


2  —   [  I  +  cos  2t,  ] 


{6e) 


E^. 


2FULD 


V 


—  (  1  +  cos  3t,  J 


(6f) 


.All  parameters  are  chosen  to  produce  the  energy  in  kcal  mol.    Drstance  is  measured  in  A 


Since  only  3A'-6  degrees  of  freedom  are  associated  with  a  molecule  of  .V  atoms,  ue  chose 


to 
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fix  6  variables  to  the  plane  r  =  0  by  means  of  an  energy  constraint  term  isofr  constramt).  This 
procedure  has  two  advantages:  it  leaves  the  problem  as  an  unconstrained  minimization  type,  and  it 
avoids  the  problem  of  an  indefinite  Hessian  due  to  translation  and  rotation  invariance.  The 
constraint  constant  shouJd  be  sufficiently  large  to  assure  that  near  a  solution  the  constraint  energy 
is  identically  zero:  with  a  good  minimization  strategy,  this  term  becomes  zero  after  only  several 
Iterative  steps.  (Note  that  one  Newton  step  is  required  to  minimize  the  quadratic  penalty  function 
alone).    Thus,  to  fix  the  first  three  atoms  to  the  v  .  y  plane    (01",  CI'   and  C2")  we  add: 

^Fix  ^  ^-^  <  ^T.i  +  -^i.2  "^  *-T.3  +  *-'2.;  "^  ^'2.3  +  ■'^5.3  >  ■  *^S' 

Several  modifications  have  been  introduced  in  our  potential  energy  function  (compare  6c. d  to 
3c).  The  bond  strain  terms  involve  a  difference  of  square'i  of  the  bond,  lengths  rather  than  a 
difference  of  bond  lengths;  the  bond  angle  strain  terms  involve  a  difference  of  cosines  of  angles 
rather  than  a  difference  of  angles.  Two  advantages  follow  from  these  modifications.  First,  the 
expressions  are  mathematically  simpler:  vector  square  roots  are  not  necessary  for  evaluating  the 
bond  strain  energy,  and  in  view  of  our  bond  angle  cosine  expression  (4).  both  Ebosd  -'"^  ^basg 
are  fourth-degree  polxnnmials  in  the  coordinate  variables.  Differentiation  of  these  term>  is  then 
straight-forward  and  the  potential  for  round-off.  a  critical  issue  for  large-scale  problems,  is 
reduced.  Second,  the  bond  angle  term  is  now  an  infinite  Taylor  series  in  powers  of  (  9  -  6  i.  For 
large  angle  deviations,  this  expression  may  happen  to  be  more  suitable  than  a  harmonic  potential. 
Furthermore,  it  is  periodic.  Indeed,  this  choice  of  bond  angle  potential  may  have  contributed  to 
the  accurate  bond  angle  values  that  were  generated  for  deoxyribose. 

For  comparison,  these  functions  can  be  matched  to  the  more  common  harmonic  strain 
functions  (3c)  by  equating  Taylor  coefficients  up  to  second  order.  The  relation  to  the  coefficients 
of  the  harmonic  potentials  (denoted  temporarily  by  .ST  and  S2')  is  then  given  by: 

S\   =  ^  S\'  7-  (  7a ) 

S2  =  S2'  sin  -  e  (7b) 
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For  a  unified  treatment  of  all  local  energy  terms  as  polynomials  of  the  cartesian  \ariahles.  ue 
express  the  torsional  potentials  (6e.t')  usmg  the  dihedral  angle  cosine  expression  (5)  in 
combination  with  the  trigonometric  identities 

cos  2t  =  2  cos  -7  -   1  .  (8a) 

cos  37  =  4  cos  '7  —  3  cos  7  .  (8b) 

W'e  now  describe  in  detail  the  choices  for  the  constants  associated  with  each  energy 
component. 


3.2a     Lennard-Jones  Interactions 

Two  different  methods  can  be  considered  for  computing  the  Lennard-Jones  parameters  A,j 
and  B,j  of  equation  (6a).  One  is  based  on  a  simple  scaling  of  the  Lennard-Jones  curve  for  each 
pairwise  interaction  [45.46]  and  the  other  involves  the  Slater-Kirkwood  equation  [47,-1-8:  see  also 
-^9-52  for  explanation  of  increased  Van  der  Waals  raddi;  and  53  for  parameters].  For  this 
calculation,  the  attractive  Lennard-Jones  coefficient  .4,,  was  obtained  from  the  .Slater-Kirkwood 
equation,  and  the  repulsive  coefficient  ff,,.  from  A,,  and  the  additional  requirement  that  the 
original  well  depth  of  the  Lennard-Jones  energy  is  maintained.  Details,  including  parameter 
tables,  are  given  in  Appendix  B. 

We  include  in  the  non-bonded  energy  component  all  interactions  between  individual  atoms 
separated  by  3  bonds  or  more.  In  cases  where  at  least  one  member  of  the  interacting  pair  is  an 
atom  :^roup.  interactions  separated  by  two  chemical  bonds  are  additionail)  considered.  Non- 
bonded  interactions  among  ring  atoms  are  not  considered  since  a  separation  of  3  chemical  bonds 
in  one  direction  involves  a  2-bond  separation  m  the  other  direction. 
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3.2b     Electrostatic  Potential 

The  eIectro^tatic  interaction  between  a  pair  of  atoms  having  charges  Q,  and  Q  and  separated 
by  a  distance  r,^  is  represented  in  this  calculation  by  the  standard  coulomb  potential  of  equation 
(6b)  modified  by  a  distance  dependent  dielectric  function  D(r,j)  =  ee  '  .  as  suggested  by 
Srinivasan  and  Olson  [54.55].  A  dielectric //(/ic-f/on  is  intr<-)duced  to  account  for  charge  shielding 
by  solvent  molecules.  For  a  neutral  molecule,  an  effective  dielectric  constant  of  e  =  4.0  has 
generally  been  used  for  small  distance  separations;  as  interatomic  distance>  increase,  this  value 
must  be  modified  to  model  penetrating  aqueous  solvent  molecules.  It  i.s  ^tlll  difficult,  however,  to 
make  a  reliable  estimate  of  the  dielectric  expres.sion.  Improved  treatments  of  solvent  effects  and 
charged  species  are  currently  being  investigated  [56.57.  for  instance]. 

Selection  of  partial  charges  has  been  a  difficult  issue  in  molecular  mechanics  calculations 
since  quantum  mechanical  approaches  yield  significantly  different  values  [58  and  references  cited 
therein].  However,  a  resolution  is  on  the  horizon  with  experimental  determination  of  these 
quantities  [59.60].  For  our  calculation,  the  partial  charges  Q,  for  individual  atoms  and  atom  groups 
are  taken  from  Olson  [17]  and  summarized  for  completeness  in  table  IV  of  .Appendix  B. 

3.2c     Bending  and  Stretching  Parameters 

Equilibrium  bond  lengths  used  are  1.52  A  for  C  —  C  bonds.  1.42  A.  for  C-0  bonds. 
1 .47  A  for  C-N  bonds,  and  1 .0  A  for  C-H.  \-H.  and  O-  H  bonds.  .All  equilibrium  bond  angles 
are  taken  to  be  tetrahedral  (  109.47"). 

In      the      choice     of    force-field      constants,      we      were      guided    by    the  general     relation 
S\  ■  L^    >  S2  >  V'2,  V'3  where  Z,  is  a  characteristic  length  iL   =    1.0    .A),     in  accordance  u  ith  the 
degree  of  flexibility  associated  with   the  geometric  quantities   -   bond  lengths,  bond  angles,  and 
dihedral  angles.    For  the  bond  length  constant  SI.  we  use  a  uniform  value  of  100.0  kcal  mol  .A 
for  all  bond  types.  For    bond  length  values  of  1.0  A  and  1.52  A  .  for  example,  the  corresponding 
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harmonic  stretching  coefficients  are  400.0  kcal  mol  A  -  and  92-1-  16  kcal  mo!  A  -.  respectively  (see 
equation  7a I.  Note  that  the  choice  oi  S\  is  not  critical,  since  any  \alue  that  can  keep  the  bond 
length  fluctuations  small  is  appropriate. 

For  the  bond  angle  constant  S2  we  use:  60.0  kcal  mol  for  all  C-O  — C.  C  — C-C  and 
O-C-C  bond  angles:  15.0  kcal  mol  for  N-C-0  .  N'-C-C  and  all  angles  involving  one  or  more 
hydrogens.  The  bending  constant  of  60.0  kcal  mol  corresponds  to  53.33  kcal  mol  deg-  in  the 
harmonic  bending  functions,  and  similarly  15.0  to  13.33  (see  equation  ^b).  These  values  v\ere 
chosen  to  be  in  the  range  of  other  parameter  values  that  have  been  used  [17.34].  By  treating  all 
bond  angles  of  a  given  type  uniformly  regardless  of  their  location  in  the  molecule,  ue  can  test 
whether  the  correct  geometry  associated  v\ith  the  puckered  structures  uill  result  from  energy 
minimization  rather  than  targeting  strategies. 

3. 2d     Torsional  Parameters  • 

Torsional  barriers  and  periodicities  are  estimated  from  experimental  data  obtained  principally 
by  spectroscopic  methods  (NMR.  IR  .  microwave.  Raman)  for  low  molecular  weight  compounds. 
.According  to  a  theory  developed  by  Pauling  [61].  potential  barriers  to  internal  rotation  arise  from 
e.xchange  interactions  of  electrons  in  adjacent  bonds  and  are  thus  similar  for  molecules  uith  :he 
same  orbital  character.  This  theory  has  allowed  tabulations  of  barrier  heights  as  class  averages 
[47,62-64].  Since  barriers  for  rotations  about  various  single  bonds  in  nucleic  acids  and  proteins 
are  not  yet  available  experimentally,  they  must  be  estimated  from  analogous  chemical  sequences 
in  low  molecular  weight  compounds. 

Twofold  and  threefold  potentials  are  used  in  the  present  calculation.  .A  threefold  torsional 
potential  exhibits  three  maxima  at  O'^  .  120'^  .  and  240'^  and  three  minima  at 
60  .  1<S0  .  and  300  .  In  ethane,  for  example,  all  maxima  are  energetically  equi\alent  and 
correspond  to  the  eclipsed  or  cis  form,  and  all  minima  correspond  to  the  energetically  equivalent 
and  preferred  staggered  or  trans  form.    In  other  molecular  sequences,  the  local  minima  ma\   not 


be  ot'  equal  energies.  The  gauche  forms  at  60"  and  300'^  may  be  higher  in  energy  than  the  trans 
form  at  180".  as  m  /!-butane  [63],  or  lower  in  energy  as  m  l-Z-ditluoroethane  and  certain 
C-C-C-O  linkages  [65].  A  combmation  of  a  twofold  and  a  threefold  potential  can  be  used  to 
reproduce  both  the  cis/trans  and  trans  gauche  energy  differences.  If  ue  denote  the  experimental 
energy  barrier  by  AV  and  the  empirical  potential  energy  function  as  given  in  (6)  by  E.  then  the 
torsional  parameters  V2  and  V'3  for  a  given  rotation  t  are  computed  from  the  following  relations; 

AV  =  F        -F 

'-^  ^  CIS  trans  '^-  =  0^         7=180" 

=    V3   +    [  Eu    +   E(^oLL    "•"   ^BOSD   ^   ^BASG  ]  7  =  0-' 

-    [  Eu   +   E(^Qi^i    +   EgQ^Q   +   ^BASC  ]  7=  ISO"  ^*^^^ 

^^  trans,  gauche  '- 7=  igQ""      t  =  AO" 

=    ~  ^-    +    [  ^U    +   ^COLL    ^   ^BO\D    +   ^'s.A.VG  1  -=180" 

-   [  Eu    -I-   Ef^oCL    '^   ^BOSD    +   ^BASG  ]  7  =  60"  *  "^^ ' 

The  formulas  above  may  appear  simple  but  in  order  to  evaluate  the  energy  as  a  function  of  t. 
cartesian  coordinates  of  the  molecule  must  be  generated  in  terms  of  t.  For  small  molecules, 
deviations  of  bond  lengths  and  bond  angles  from  equilibrium  values  as  a  function  of  x  can  be 
obtained  experimentally,  but  when  torsional  parameters  are  computed  for  large  molecules,  these 
deviations  are  difficult  to  approximate.  .A  Mmpiificatu)n  can  be  made  by  assuming  fixed  bond 
lengths  and  bond  angles  and  by  calculating  only  non-bonded  energy  differences.  In  theory,  then, 
every  different  parametrization  of  the  non-bonded  coefficients  requires  an  estimation  of  the 
torsional  potentials  V2  and  V'3  to  produce  a  consistent  set. 


For  the  present  investigation  of  deoxyribose.  we  follow  Olson  in  assigning  the  values  used 
for  V'3:  2.8  kcal  mol  for  rotations  about  C-C  bonds  and  1.8  kcal  mol  for  rotations  about  C-O 
bonds  [17j.  These  values  were  obtained  to  reproduce  the  experimentally-observed  barriers  for  n- 
butane  (3.5  kcal/mol)  and  inethyl  ethyl  ether  (2.53  kcal,  mol).  respectively,  together  s\  ith  the 
Lennard-Jones  potential.  The  values  used  for  V2.  namely  0.2  kcalmol  for  all  O  — C-C-C 
rotations  and  1.0  kcalmol  for  all  O-C-C-O  rotations,  were  introduced  to  reproduce  the 
trans  gauche  energy  differences  for  similar  molecular  sequences  [17]. 


.Additivity  of  torsional  energies.  Torsional  potentials  discussed  so  far  are  determined  to 
reproduce  rotational  barriers  for  only  one  rotation  for  a  given  4-atom  sequence  in  a  gi\"en 
molecule.  .As  was  first  mentioned  by  Harvey  and  Prabhakaran  [191.  previous  treatments  of 
dihedral  rotations  were  somewhat  arbitrary  since  only  a  subset  of  all  rotations  were  chosen  to  be 
modeled  in  the  energy  function.  When  a  dihedral  angle  about  ti  B-C  bond  can  be  defined  by 
more  than  one  .A  — B-C  — D  sequence,  several  issues  must.be  resolved: 

1 )  Treatment  of  threefold  potentials. 

How  should  the  V'3  parameter  be  selected  and  adiusted'.' 

2)  Treatment  of  twofold  potentials. 

These  potentials,  as  used  in  this  study,  are  different  from  the  threefold  terms  because  their 
force  constants  V2  are  associated  with  a  4-atom  sequence  and  not  a  2-atom  bond.  Thus, 
unlike  the  threefold  potentials,  the  chemical  type  of  atoms  .A  and  D  in  the  ,A-B-C-D 
sequence  affects  the  choice  of  torsion  parameter.  When  more  than  one  such  sequence  is 
associated  with  the  same  atom  bond,  should  the  V2  parameter  be  adjusted? 

3)  Combined  twofold  and  threefold  potentials. 

Fs  rurther  adjustment  required  when  the  same  dihedral  ansle  is  included  in  b(nh  terms'.' 


Since  torsional  potentials  are  semi-empincal  in  nature,  the  answers  to  these  questions  are  not 
completely  clear.  By  considering  equations  (9a)  and  (9b).  the  answer  to  question  3  may  appear 
simple:  no  conflict  arises  from  the  inclusion  of  one  particular  rotation  in  both  the  twofold  and 
threefold  potentials,  since  the  two  parameters  V'2  and  r3  are  determined  by  two  separate 
equations.  However,  computing  the  cis  trans  and  trans-gauche  energy  differences  for  dijfereni 
molecules  is  not  ideal.  For  example,  in  the  specific  parameter  values  described  above,  rotations 
about  the  C-C  bond  were  chosen  to  reproduce  the  cis  trans  barrier  for  /!-butane  and  the 
trans  gauche  barrier  for  various  C-C-C  — O  and  0-C-C  — O  sequences.  Nonetheless,  this 
procedure  may  be  justifiable  since  we  are  trying  to  incorporate  two  separate  effects  observed  in 
small  molecules  in  the  study  of  larger  and  more  complicated  systems. 

Question  1  was  already  raised  in  the  molecular  dynamics  study  of  ribose  by  Harvey  and 
Prabhakaran  [19].  In  their  study,  a  compromise  modeling  between  the  minimum  set  of  5 
endocyclic  dihedral  angles  and  the  ma.ximum  set  of  all  16  heavy  atom  rdtations  was  implemented. 
For  our  investigation  of  deoxyribose.  we  examine  the  modeling  of  all  12  dihedral  rotations 
associated  v\ith  heavy  atoms  or  atom  groups.  To  reproduce  a  particular  cis  trans  energy 
difference  about  a  particular  bond,  1'3  is  adjusted  by  a  simple  division  by  the  number  of  modeled 
rotations.  Our  results,  described  in  the  next  section,  indicate  that  this  simple  modification  of  the 
torsional  parameters  can  allow  a  systematic  treatment  of  all  rotations  in  a  given  molecule.  W'e  list 
all  the  rotations  and  their  as.sociated  torsional  parameters  in  Table  I. 

Question  2  appears  more  difficult.  For  illustration,  consider  all  4  rotations  about  the 
C3'  — C4'  bond  as  listed  in  Table  I.  Four  threefold  rotations  are  associated  about  this  endocyclic 
bond  (with  V'3  of  2.8  kcaLmol  equally  distributed  among  them),  as  are  3  twofold  rotations:  2  of 
which  are  C-C-C-0  sequences  and  one  of  vshich  is  0-C  — C  — O.  Since  the  4-atom  sequences 
involve  different  atoms,  it  seems  that  no  V2  adjustment  is  necessary.  However,  this  conclusion 
appears  inconsistent  when  we  consider  the  above  apportionment  dcMsed  for  obtaining  the  V'3 
parameters  for  the  threefold  potentials.  We  reserve  this  question  for  future  thought  and 
experimentation.     In  our  present  study,  the   gauche  contributions  are  quite  small  and  therefore 


each  4-atorn  sequence  can  be  individually  considered  at  full  weight. 


33      Organization  for  Minimization 

The  form  of  the  potential  energy  function  uas  constructed  to  allow  direct  differentiation  uith 
respect  to  the  coordinate  variables  and  to  ensure  that  the  form  of  the  derivatives  is  as  Mmpie  as 
possible.  Organization  of  the  data  Ntructures  that  store  the  geometric  information  is  also  crucial  to 
the  efficiency  of  the  minimization  strategy.  The  bond  lengths,  bond  angles,  dihedral  angles  and 
associated  parameters  are  stored  in  arrays  that  facilitate  accessing  the  data  for  the  function  and 
derivative  evaluations. 

Three  stages  comprise  the  con.struction  of  such  data  tables.  First,  the  individual  atoms  are 
numbered  and  are  assigned  identification  labels  for  atom  types;  then  the  individual  bonds  are 
numbered.  Second,  two  connectivity  arrays  are  entered  to  specify  bonded  atoms  and  neigbors  of 
each  atom,  and  four  parameter  arrays  to  specify  stretching,  bending  and  torsional  constants 
associated  with  general  chemical  sequences.  Third,  a  subroutine  is  called  to  construct  lists  of  all 
the  bonds,  bond  angles  and  dihedral  angles  in  the  molecule  and  then  to  associate  the  proper 
energy  parameters  to  each  sequence.  In  this  way.  changing  energy  parameters  is  an  easy  task,  and 
calculation  of  the  energy  terms  and  derivatives  requires  a  simple  reading  "down  the  list"  where  all 
the  data  have  been  preprocessed. 

The  required  derivatives  are  assembled  in  a  modular  fashion.  Function  subprograms  are 
used  to  compute  the  derivatives  of  various  geometric  quantities  with  respect  to  the  coordinate 
variables    r^  /  .    for    k  =  ] V     and     /  =   1  .  2  .  3  .    For    example,    the    deri\ative    of    the 

interatomic  distance  vector  r,^  with  respect  to  .v;.  ,.    - — -.    is  siven  bv  the   function  deni.j .k .1 ) . 

2n*J     :i        -,  hy  dderi  i .)  .k .1  .m  .n\ .    The  expression   for  a  vector  inner  product  insolvine  4 

atoms    (X,  -  \.)  ■  (X,-  -  x^)      is    given    in    function    proddj  J'  j' )    and    the    derivatives    of    this 

d  prndt  i.j  j'  .j'  ]                         d'  prndii  j  i'  ;') 
quantity.  -  and       ' '-^ .        m       the       function       subnroerams 


depii  .j  A' .]'  .k.l)  and  ddepii  .j  A'  .j'  .k  J  .in  m)  .  respectively.  Similarly,  a  function  cnsine(  i ./  A'  J') 
evaluates  the  dihedral  angle  cosine  approximation  given  in  equation  (5i.  and  its  associated 
derivatives  are  given  in  two  subprograms. 

This  construction  allows  products  and  quotient.s  in  the  energy  expression  to  be  differentiated 
simply  and  elegantly  and  can  reduce  the  potential  for  errors.  With  analytic  first  and  second- 
derivatives  at  hand,  Newton  methods  can  be  implemented.  The  two  particular  Newton  algorithms 
that  we  have  tested  are:  Gill  and  Murray's  modified  Newton  method,  and  a  truncated  Newton 
method  adapted  for  potential  energy  minimization.  The  former  sohes  the  Newton  equation 
directly  by  the  modified  Cholesky  factorization  and  is  thus  suitable  for  functions  of  about  100 
variables  or  less.  The  latter  solves  the  Newton  equation  iteratively  by  a  truncated  preconditioned 
CG  and  is  thus  attractive  for  large-scale  problems.  Full  details  of  the  algorithms  are  pro\ided  in 
the  accompanying  paper  [68]. 

All  computations  described  in  the  next  section  were  performed  on  4  Vax  S600  at  Courant 
Institute.  New  York  Universitv. 


4.     RESULTS  AND  DISCUSSIONS 

The  modeling  and  minimization  techniques  developed  m  this  study  were  designed  tor 
unconstrained  potential  energy  minimization  in  cartesian  coordinate  space.  Minimized  structures. 
however,  are  only  as  good  as  the  chosen  parameter  set.  Furthermore,  the  conformation  space  uill 
be  surveyed  properly  only  with  powerful  minimization  strategies.  Thus,  to  evaluate  the  reliability 
of  the  results,  we  must  ensure  that  only  two  \oca\  minima  are  obtained  for  the  deoxyribose  model 
and  that  their  geometries  and  energies  are  consistent  with  observed  data.  For  nucleic  acid  sugars, 
even  these  simple  criteria  have  provided  a  continuing  challenge.  Going  a  step  beyond  consistency. 
v\e  were  interested  m  examining  the  effects  of  changes  in  the  energy  function  form  and 
parameters  on  the  structures  obtained  and  on  the  energy  differences  among  them.  In  particular. 
the  effect  of  the  gauche  potential  and  different  choices  for  modeled  dihedral  angle  sets  modeled  in 
the  threefold  torsion  potential  are  important  for  systematic  parameter  assignments  in  DN.A 
models. 

.A  convenient  way  to  study  all  these  issues  is  to  formulate  several  parameter  sets  and  examine 
results  from  energy  minimization  and.  additionallv .  from  unminimized  Energy  vs.  P  (the 
pseudorotation  parameter)  curves.  In  previous  pseudorotation  investigations.  E  vs  P  profiles 
were  usually  generated  using  artificial  constraint  terms  in  combination  with  minimization 
[18.20.34].  These  constraints  restricted  the  backbone  dihedral  angle  passing  through  the  ring,  all 
5  endocyclic  dihedral  angles,  or  one  endocyclic  dihderal  angle,  to  assume  fixed  values  associated 
with  the  pseudorotation  path.  Clearly,  results  are  constraint  dependent;  the  stronger  the 
constraint,  the  higher  the  energy  barriers  expected.  Since  full  cartesian  minimi/jtion  is  our  major 
computational  tool,  results  from  E  vs.  P  curves  that  are  unminimizcd  are  sufficient  lor  examining 
influences  of  different  parameter  sets  and  selecting  the  best  one.  Thus,  the  E  \sP  prntiles  are  not 
intended  to  determine  barrier  heights  as  these  will  be  exaggerated  uithout  complete  cnerg\ 
relaxation. 

The  minimization  results  also  provide  an  opportunity   to  examine  the  deviation  from  "ideal" 
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pseudorotation  behavior.  We  will  analyze  the  approximation  made  in  the  pseudorotation  equation 
(2)  for  the  5  endocyclic  dihedral  angles  by  two  different  comparisons.  First,  to  compare  P  and 
"ma.x  corresponding  to  our  minimized  structures  uith  the  values- determined  for  analyzed  crystal 
structures,  we  compute  P  and  T^a.^  from  an  exact  Fourier  series  (described  in  the  next  section). 
Second,  to  examine  the  analytic  fitting  provided  by  (2)  for  the  individual  dihedral  angles,  we 
compute  the  sum  of  squares  of  the  angle  deviations  between  the  values  generated  by  minimization 
and  the  values  predicted  by  (2)  u  ith  the  computed  Fourier  Series  parameters  P  and  T^n^x  ^^  c  use 
the  Fourier  series  representation  because  it  provides  an  equal  treatment  of  all  dihedral  angles.  .A 
simpler  approximation  for  the  pseudorotation  parameters  can  be  obtained  from  the  formula 
derived  by  .Altona.  Geise  and  Romers  [291  from  the  pseudorotation  equation  (2): 


T4  +  Ti-T,-T,| 

tan  P  =  ^ 


2  Tt  (sin  777  +  sin  ^j 

10  D 


■After  computing  P  from  this  formula.  T^ax  can  be  obtained  from  (2)  for  j  =  2.  The  disadvantage 
of  this  procedure  is  that  it  is  somewhat  arbitrary:  a  different  dihedral  angle  labeling  scheme  uill 
produce    different  p.seudorotation  phase  angle  and  phase  shift. 

The  present  section  is  organized  as  follows.  First,  we  will  describe  the  different  energy  sets 
examined,  the  details  of  generating  coordinates  as  a  function  of  P.  and  the  expression  of  the 
Fourier  series  and  derived  pseudorotation  parameters.  Second,  we  will  describe  the  structural, 
energetic  and  geometric  details  for  the  deoxyribose  conformations  generated  by  minimization,  and 
then,  for  these  structrues.  ue  will  analyze  the  dihedral  angles  in  the  context  oi  pseudorotation. 
Third,  we  will  briefly  analyze  the  obtained  E  vs.  P  profiles.  Fourth,  ue  uill  conclude  with  a 
comparison  of  the  tuo  approaches  —  energy  minimization  and  E  vs.  P  curves. 
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4.1     The  Potential  Energy  Functions  Investigated  and  Methods  of  Analysis 

Eight  different  variations  of  the  potential  energy  function  described  in  Section  3  are  examined 
(see  Tables  I  and  II).  The  non-bonded  parameters  remain  constant  throughout  the  calculations. 
but  different  choices  and  combinations  of  bond  length  constants,  bond  angle  constants,  and 
torsional  terms  comprise  the  various  energy  sets. 

The  first  potential  energy  function  Ntudied  is  composed  of  non-bonded  interactions,  bond 
length  and  bond  angle  terms,  and  a  threefold  torsional  potential  for  the  5  endocvclic  dihedral 
angles  only.  For  the  choice  of  bond  angle  bending  constants,  we  wanted  to  ensure  that  the  correct 
puckering  geometrv  emerges.  Therefore,  it  uas  not  clear  whether  the  tetrahedral  angle  tormed  b>' 
a  3-atom  sequence  mvolving  at  least  one  hydrogen,  should  be  made  to  stay  \ery  near  to 
tetrahedral  or  be  made  extremely  flexible  to  accomodate  the  endocyclic  bond  angle  arrangements. 
Thus,  in  set  1  we  tried  a  verv  larae  bendins  constant  for  these  anales  -  200.0    kcal  mol.  and  in  set 

-  -  T 

2.  a  much  smaller  value  of  15.0  kcafmol.    Since  results  were  qualitatively  similar,  we  adopted  the 
latter,  more  intuitive  choice  in  the  remaining  energy  parameter  sets  (2-8). 

We  then  added  a  twofold  or  gauche  potential  to  the  energy  function  (set  3).  a  weighted 
threefold  torsional  term  to  include  all  heavy  atom  rotations  about  C-C  and  C-0  bonds  (set  4), 
and  then  examined  their  combined  affect  (set  5).  In  sets  6-8.  we  tried  various  combinations  of 
bond  length  constants,  bond  angle  constants,  and  torsional  terms  in  order  to  examine  the  overall 
relationship  between  the  structures  obtained  and  the  potential  energy  form. 

.Minimization  runs  were  started  from  points  corresponding  to  P  values  at  18'  intervals  along 
the  pseudorotation  cvcle  and  at  randomly  perturbed  points  from  those.  Cocudinates  as  a  function 
of  P  were  generated  using  Pearlman  and  Kim's  procedure  [66].  for  the  skeletal  ring  atoms. 
Exocyclic  substituents  were  attached  to  ticcupy  positions  as  close  to  tetrahedral  in  a  local  bisecting 
ccmrdinate  system  (see  Figure  5i.  Note  that  this  procedure  was  only  used  to  obtain  the  initial 
guess  to  start  the  minimization.  Once  started,  the  minimization  procedure  makes  no  use  of  the 
pseudoroation  cycle.    Results  for  minimization  runs  are  reported  in  Table  III. 
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For  the  E  vs.  P  profiles,  the  coordinate  generation  procedure  mentioned  above  uas  used  at  1 
inter\al>.  and  the  variation  of  the  total  energy  and  individual  energy  components  were  plotted  as  a 
function  of  P.    The  curves  for  6  representative  energy  sets  are  shown  in  Figure  7. 

For  the  minimization-generated  structures,  analysis  of  the  pseudorotation  approximation  is 
accomplished  by  considering  the  exact  representation  of  the  5  endocyclic  dihedral  angles 
T,  .  /  =  0 4.  given  bv: 


4  I  ^^  j  k 

T^  =    V   a;,  f     '  ;■  =  0 4  .  (  lOai 

k  =  \) 

The  5  Fourier  coefficients  are  2iven  bv  the  inverse  Fourier  series 


1        4  -i^jk 

Oj  =  -  ^   T^e       '  /t  =  0  ,  .  .  .  ,  4  .  (10b) 

^    7  =  0 

Representing  each  complex  Fourier  coefficient  by  a^.  =  A^.  e  "  .  we  obtain  by  periodicity 
a,--;,  =  a^  (the  bar  denotes  the  complex  conjugate),  and  thus  the  series  can  be  simplified  to  the 
cosine  expansion 


T,  =  An  +  2,4,  cos  (  9,  +  ^^  )  +  2A^  cos  (  9.  +  — r^  )  .  (11) 

where 


A„  =  7  2   Tj  .  (11a) 


"*     '  l/lT/ 


,4,  cos9,  =    —  S   '  T,  COS— T-^  )  .  /  =   1  .2    .  (lib) 

-     ;=0 

1      -"  li-ni 

A,  sin9,  =  -—  y   (7,  sin — r^  )  .  /  =  1  .2    .  (1  Ic) 

By  comparing  equation  (11)  with  (2).  it  is  clear  that  the  pseudorotation  approximation  coincides 
with  the  analytic  Fourier  >eries  representation  when  Ay  =  A^  =  0.    P  and  x^^^  correspond  to  the 
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two  Fourier-derived  expressions: 


P     ~     Q.  +  ^  'iri'i  "max     ~     --^:-  (1^) 


4.2     Results  from  Minimization 

For  ail  energy  parameter  sets  and  for  all  initial  molecular  conformations,  only  tuo  local 
minima  are  obtained.  As  will  be  discussed  below,  the  minima  correspond  closely  to  the  ideal 
C2"-endo  and  C3"  — endo  sugar  puckers.  The  North  pucker  is  obtained  from  >tarting  structures 
corresponding  to  increasing  P  values  of  2''0"  to  72".  and  the  South  pucker  is  obtained  from  a  P 
range  of  90'^  to  252  (see  Figure  bi.  Random  starting  points  predict  one  of  the  two  identical 
structures.  These  results,  as  we  see  most  clearly  from  Figure  6.  support  the  idea  that  the 
pseudorotation  surface  is  a  low-energy  path. 

Quadratic  convergence  is  realized  near  a  local  minimum  for  both  the  modified  Neuton 
method  and  the  truncated  Newton  method.  For  both  Newton  methods,  about  IS  Newton 
iterations  and  21  function  and  gradient  evaluations  were  required  to  reach  the  convergence 
criterion  of  a  gradient  norm  of  less  than  1.0  x  10" '.  .Additionally,  each  Neuton  iteration  required 
one  evaluation  of  the  analytic  Hessian  for  the  modified  Newton  method  and  one  Hessian  \ector 
multiplication  in  the  truncated  Newton  method.  Thus,  for  an  optimization  problem  of  this  size, 
performance  of  Newton  methods  based  on  a  direct  solution  a\  the  Newton  equation  is  very 
similar  in  computational  effort  and  convergence  properties  to  performance  o(  Newton  methods 
based  on  an  iterative  linear  solver.  The  truncated  Newton  will  undoubtedly  be  superior  for 
functions  with  a  larger  number  of  variables  when  the  search  m  the  large  conformational  space  will 
be  more  adaptive  to  the  progress  made  and  more  directed  toward  the  important  conformational 
regions.  .Moreover,  the  storage  requirements  of  the  truncated  Newton  method  can  be  made  less 
than   the   modified   Newton   method   for  any   problem   size   by   using  a   preconditioned  Conjugate 


Gradient  variant  that  requires  analytic  evaluation  and  storage  of  only  the  "local"  Hessian  elements 
entering  from  the  bond  lengths,  bond  angles  and  torsional  potentials.  In  other  words,  the  long- 
range  non-bonded  terms  in  the  Hessian  matrix  need  not  be  computed  at  all!  (details  are  provided 
in  the  accompanying  paper  [68]).  .Application  of  the  truncated  Newton  method  to  large  potential 
energy  functions  is  promising,  because  our  experience  with  the  truncated  Newton  code  on  various 
test  problems  suggests  that  function  form  rather  than  problem  size  governs  the  algorithm'^ 
performance.  We  are  currentlv  working  on  efficient  implementation  of  the  method  to  larger 
molecules. 

STRUCTURES.  We  summarize  in  Table  III  the  generated  endocyclic  bond  angles,  selected 
dihedral  angles,  and  the  values  of  all  energy  components  for  the  two  structures  obtained  from  all 
minimization  trials.  For  comparison,  we  list  corresponding  bond  angles'  for  the  C2"  — endo  and 
C3"-endo  conformations  obtained  experimentally  [7.9]  and  endocyclic  dihedral  angles  predicted 
theoretically  by  the  pseudorotation  equation  (2)  with  "ma.x  ~  ^^^"  •  Additionally,  the 
pseudorotation  parameter  P  and  puckering  amplitude  ~j^-^y_  are  computed  for  each  structure  from 
the  coefficients  of  a  Fourier  series  by  equation  (12l.  These  values  are  used  to  compute  the  5 
dihedral  angles  from  (2)  and  are  compared  to  the  dihedral  angle  value  generated  by  minimization 
by  summing  the  squares  of  the  individual  angle  deviations.  This  value  is  an  indication  of  the 
approximation  involved  in  pseudorotation. 

The  estimated  pseudorotation  angles  for  the  minimization-generated  structures  range  from 
9'  -  25'  with  an  average  (over  all  energy  parameter  sets  I  of  19''  for  the  North  pucker 
(C3'-endo);  the  estimated  P  values  range  from  158*"'  -  16S''  and  average  16l"  for  the  South 
(C2"-endo).  These  values  are  very  close  to  the  ideal  phase  angles  of  18"  and  162".  The 
puckering  amplitude  Tj^^^  has  an  average  c^f  40.l"  for  the  North  and  38.5"  for  the  South  puckers. 
This  is  clearly  a  small  difference,  but  such  a  trend  is  observed  consistentK  for  all  different 
potential  energy  sets  and.  moreover,  has  been  noted  previously  in  theoretical  investigations 
[19.34]   and   experimental  analyses    [6,7].   Thus,   conformational    interchanges   by   pseudorotation 


with  variable  amplitudes  of  puckering,  or  through  a  planar  contormation.  may  occur. 


ENERGIES.  The  small  but  significant  energy  differences  between  the  North  and  South 
conformations  illustrate  the  dilemmas  that  arose  with  regard  to  reproducing  the  South  preference 
for  deoxynbose  fragments  in  earlier  potential  energy  investigations  of  deoxyribose.  ClearK.  the 
only  energy  functions  that  predict  a  lower  energy  state  at  C2'-endo  are  those  that  include  the 
gauche  potential  ( >ee  Tables  III  and  IV). 

The  non-bonded  and  threefold  torsional  energy  components  are  re>pecti\ely  about  0.4 
kcafmol  and  0.3  kcalmol  lower  for  the  North  pucker  than  for  the  Si)uth  pucker,  while  the  bond 
length  strain  terms  are  small  and  nearly  equal.  In  combination  with  the  bond  angle  component 
v\hich  IS  about  0.4  kcafmol  lower  for  the  South  conformer.  the  total  energ_\  is  lower  for 
C3"-endo.  Thus,  only  the  gauche  potential  can  sufficiently  penalize  the  energy  at  the  North 
minimum  to  make  the  South  or  C2'-endo  a  global  minimum. 

The  gauche  term  is  dominated  by  the  0 1"  — C4'-C3'-03  rotation.  The  endocyclic 
O  — C-C-C  sequences  contribute  about  the  same  energy  to  both  local  minima,  but  the 
O  1  ■-C4"-C3"-03"  sequence  exhibits  a  major  conformational  difference.  This  dihedral  angle  is 
gauche  I  ~  96")  at  the  South  pucker  and  trans  (  ~  157")  at  the  .North.  Since  the  barrier  height  is 
higher  for  0-C-C-O  than  for  O  — C-C-C  rotations,  the  gauche  contribution  from  this 
sequence  alone  lowers  the  South  pucker  energy  by  about  0.8  kcafmol.  The  C  r-C2'-C3'-03" 
sequence  is  in  the  gauche  range  in  the  South  and  the  trans  range  in  the  North,  but  its  effect  is 
balanced  by  the  reverse  trend  for  the  C5"-C4"-C3"-03"  sequence.  The  observed  correlation 
between  0 1  •-C4--C3--03"  (~  t,  -  120")  and  C5--C4--C3"-03"  [^  -;,  +  120"i  rotation 
sequences  and  the  ring  pucker  mode  has  been  used  to  characterize  the  state  of  the  sugar  pucker  in 
the  nucleotide  backbone  [34]. 

Other  variations  in  energy  potentials  have  little  effect  on  the  relative  energy  differences 
between  the  North  and  South  minima.    L'sing  a  threefold  torsional  potential  for  all  hea\\  rotation 


sequences,  as  described  in  Table  I.  slightly  lowers  the  energy  difference  between  the  North  and 
South  conformers.  From  a  careful  examination  of  the  values  of  all  modeled  dihedral  angles,  it 
becomes  apparent  that  the  contribution  from  the  endocyclic  angles  dominates  over  that  of  other 
angles  in  the  group  with  v\hich  the  1'3  parameter  is  distributed.  Consequently,  the  overall  effect 
of  including  all  heavy  rotations  of  the  molecule  in  a  weighted  threefold  torsional  potential  is  to 
lower  the  torsional  energy  component.  As  we  will  see  from  the  Energy  versus  P  curves,  the 
threefold  torsional  potential  is  lowered  consistently  for  the  entire  pseudorotation  cycle,  and  thus 
does  not  change  the  relative  energies  between  different  sugar  puckering  modes. 

When  changing  stretching  and  bending  parameters,  ue  also  observe  very  small  relative 
energy  differences  between  the  two  minima.  Reducing  the  bond  stretching  constant  .S'l  from  100.0 
to  25.0  kcal  mol  X  "  (compare  sets  4  and  6i  changed  the  energy  difference  by  only  0.03  kcal  mol. 
Using  a  constant  bending  constant  of  60.0  kcalmol  for  all  bond  angles  of  the  molecule  (compare 
set.s  6  and  7)  changed  the  energy  difference  by  only  0.01  kcaf'mol. 'Sirnilarly.  changing  both  .S'l 
and  52  significantly  (compare  set.s  5  and  8)  changed  the  energy  difference  by  only  0.01  kcafmol. 
These  observations  suggest  that  the  model  is  qualitativelx  good  -  the  bending  and  stretching 
constants  are  appropriately  chosen  in  relation  to  each  other  and  to  other  parameters  of  the 
potential  energy  function  -  so  that  the  non-bonded  interactions  are  allowed  to  govern  the  ring 
geometry  even  with  significant  parameter  changes  in  other  energy  compiments. 

Thus,  parameter  sets  5  or  8.  both  which  include  a  gauche  potential  and  a  torsional  potential 
for  all  heavy  rotations  in  the  molecule,  are  successful  for  reproducing  the  proper 
C3"-endo/C2"-endo  relative  energy  difference. 

PUCKERED  GEOMETRY.  The  geometric  quantities  obtained  from  energy  minimization 
are  in  very  good  agreement  with  available  experimental  data  for  all  energy  sets  tested:  the 
endocyclic  bimd  angles  are  mostly  uithin  I  of  observed  values.  It  should  be  noted  however,  that 
geometric  data  available  for  deoxyribose  come  only  from  crystal  structures  o(  nucleosides  and 
nucleotides.  In  isolated  form,  the  2"  — deoxvribose  molecule  is  unstable    —  the  furanose  rina  tends 
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to  adopt  a  pyranose    or    6-membered  ring  form  [6' 


The  t'ive-membered  ring  puckers  from  a  planar  conformation  in  order  to  minirni/e  the  non- 
bonded  energy.  The  non-bonded  energy  is  greatest  in  the  planar  ring,  uhen  all  endocyclic 
substituents  are  eclipsed.  Any  puckered  form  can  reduce  this  strain,  but  the  best  arrangement  is 
obtained  when  C2"  or  C3"  puckers  out  of  the  plane  of  the  other  4  atoms  on  the  same  side  of  C5'. 
No  constituents  remain  eclipsed  at  that  position.  Since  rotational  barriers  about  C  — O  bonds 
are  lower  than  those  about  C-C  bonds,  the  torsional  energy  is  lower  uhen  endocyclic  dihedral 
angles  about  C-0  bonds  are  O'  (atoms  are  coplanar)  rather  than  about  C-C.  Thus,  the  two 
conformations  associated  with  t,)  =  O"  and  -4  =  O"  with  all  other  rotations  (Tj  .t^  and  t;) 
maximally  eclipsed  are  exacth  the  C3"  — endo  and  C2"  — endo  sugar  puckering  forms. 

As  the  ring  deviates  from  planarity.  endocyclic  bond  angles  must  necessarily  deviate  from 
planar  ring  angles  of  108'  and  from  tetrahedral  angles  of  109.47'-  which  are  inconsi.stent  with 
ring  closure  [10].  The  C-0  — C  bond  angle  remains  very  close  to  tetrahedral  since  the  ring 
oxygen  is  farthest  from  the  puckered  C2"  or  C3'  atoms  and.  furthermore,  has  no  exocNclic 
subtituents.  Bond  angles  at  the  neighboring  CI"  and  C4"  atoms  are  reduced  from  the  tetrahedral 
value  by  about  2"  to  5".  while  the  bond  angles  at  C2"  and  C3'  are  reduced  the  most  -  by  ft"  to  S'''. 
The  bond  angle  at  the  puckered  atom  (C2"  or  C3")  is  always  the  lowest  among  all  endocyclic 
angles. 

All  these  trends  are  evident  in  the  energy-minimized  structures.  It  is  interesting  that  the 
observed  trend  in  the  endocyclic  bond  angle  values 

c-o-c   >  C-C-0    >  c-c-c 

is  accurately  reproduced  as  a  result  of  energy  minimisation  and  not  as  a  consequence  of  different 
bending  force  constants  for  t.iese  angles:  all  endocyclic  bond  angles  were  assigned  a  bonding 
constant  of  60.0  kcal/mol.s  in  all  energy  parameter  sets. 

With  regard  to  bond  lengths,  the  large  bond  stretching  constant  .S"I  enforces  values  \ery  close 
to    the    observed    equilibrium    bond    lengths.    Nevertheless,    in    most    computed    geometries    the 


0r-C4'  bond  is  slightly  shorter  than  the  01'  -  CI'  bond  (  by  about  0.004  X  ).  clearly  smaller  but 
indicative  of  the  trend  observed  in  nucleoside  crystal  structures  [17].  Bond  C2"-C3'  is  also 
slightly  shorter  than  other  C  —  C  bonds,  and  computed  C  — H  bond  lengths  are  usually 
somev\hat  stretched  from  1.0  .■^. 

From  the  sum  of  squares  of  deviations  in  the  generated  endocyclic  dihedral  angles  and  those 
obtained  from  (2)  with  the  calculated  P  and  t,^^,^  (see  last  row  of  Table  III),  we  see  a  range  from 
0.7  to  60.0.  For  all  but  2  of  the  16  structures  analyzed,  the  sum  is  not  large.  The  sum.  howe\er. 
is  significant  for  the  North  minimum  of  set  2  and  the  South  minimum  of  set  6:  36.0  and  60.0. 
respectively.  In  these  cases.  Aq  and  .4i.  the  magnitudes  of  both  neglected  Fourier  coefficients  uq 
and  ill.  are  large.  Thus,  a  good  fit  of  the  pseudorotation  approximation  relies  on  ,4,i  and  ,4  i  being 
of  negligible  values.  From  our  analysis,  these  magnitudes  are  considerably  smaller  in  relation  to 
that  of  Ai  but  are  not  negligible.  Thus,  although  in  general  the  pseudorotation  description 
produces  a  good  approximation  for  P  and  x^ax-  it  does  not  always  produce  a  good  approximation  ■ 
to  the  individual  dihedral  angles. 

Mathematically,  an  explanation  for  the  small  magnitudes  of  the  neglected  lower  order  Fourier 
coefficients  can  be  provided  by  considering  the  original  Cremer  and  Pople  formulation  [Ui.  In  the        j 
construction  of  the  Fourier  coefficients  for  a  given  puckered  form,  atomic  displacements  of  the  5 
ring  atoms  from  a  chosen  reference  plane  are  given  exactly  by  a  truncated  Fourier  series  (I),  since 
this  plane  is  chosen  to  guarantee  that  A^  =  ,4;  =  0. 

In  view  of  equation  set  (  I  1 1 .  Cremer  and  Pople  have  oriented  the  plane  r  =  0  so  that  the'  ring 

coordinates    x,    are    translated    to    x,    (/  =  0 4).    whose    r-coordinates    f.v,  ■;)    satisfv    the 

J      J  J  ..^ '  . 

conditions: 

S  -^.Z  =  0  .  I   ^j.2  cos^  =  0  .  J^    v^,;  sin^  -  0  .  (  13) 

Physically,  these  conditions  imply  that  the  origin  of  the  mean  plane  is  the  center  of  mass.  and.  in 
the    case    of   small    puckering    amplitude    for    a    regular    pentagon,    that    angular    momentum    is 
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conserved  [14]. 

The  new  coordinates  x,  are  defined  relative  to  a  new   origin  so  that  (  1  la)     is  Natistied: 

x;  =  X,  -  7   ix,  .  ;•  =  0 4.  (14) 

-    k  =  0 

and  the   coordinate    r-a\is    is    defined    to    be    in    the    direction  of  the    unit    normal    n: 

.v^'.3  =  x,  ■  n  ,  ;  =  0 4  .  (15) 

uhere 

tlxt2  ,         J,         .     2tt;  ^  ,        ^  --^i 

n  =  Trz Tn   •  tl  =    >  \,  Mn — r-  .  and  t2  =    >  x,  cos — r^ 

I|tl^t2|i  -^^  ^  -,^ 

By  construction,  tl  and  t2  are  perpendicular  to  n  .  Thus,  the  inner  products  ot  tl  and  t2 
with  n  guarantee  that  the  second  and  third  conditions  of  (13)  hold  also.  This  derivation  can 
explain  s^hy  the  truncated  Fourier  series  (2)  provides  a  reasonable  approximation:  Dunitz 
has  shown  that  (2)  follows  from  (  1)  for  a  regular  pentagon  and  wuh  infinitesmal  deviations  from  a 
planar  conformation  [10]. 

Differences  between  the  Cremer  and  Pople  description  and  the  Altona  and  Sundaralingam 
description  will  occur  for  the  furanose  ring.  >ince  Dunitz'  derivation  is  not  valid  for  finite 
puckering  amplitudes  and  irregular  pentagons.  Furthermore,  the  rather  ingenious  mean  plane 
construction  of  Cremer  and  Pople  has  a  subtle  flaw:  this  plane  does  not  remain  constant  as  the 
ring  follows  the  pseudorotation  cycle.  Consequently,  even  for  a  regular  pentagon.  an\  deviation.s 
from  an  initially-chosen  plane  necessarily  lead  to  approximations  in  equation  (\),  and  in  turn,  to 
approximations  in  (2). 


4.3    Results  from    Energy    vs.     P    Curves 


STRUCTURES.    The  two  local  minima  obtained  in  the  Energy  v>.  P  curves  for  all  parameter 
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sets  occur  at  an  average  value  of  P  =  14.'-'  for  the  North  and  P  =  173'^  for  the  South,  not  as  close 
to  the  ideal  values  of  IS"  and  162''  as  those  predicted  by  energy  minimization.  .A  11  East  maxima  - 
corresponding  to  the  Ol'-endo  pucker  -  occur  near  P  =  90''.  and  all  West  ma.xima  -  the  Oi"-e\o 
pucker  -  occur  near  P  =  276'\  The  non-bonded  energies  are  high  at  these  confomations  as  a  result 
of  the  eclipsed  C2"  and  C3"  substituents  in  the  East  and  the  equatorial  orientation  of  the  base  with 
the  exocyclic  C5"  substituents  in  the  West.  Since  conformations  were  rigid  for  a  given  value  of  P  . 
barrier  heights  are  exaggerated.  Despite  this  procedure,  the  average  barrier  for  interconversion 
between  the  C2"-endo  and  C3"-endo  conformations  through  the  East  region  of  the 
pseudorotation  path  -  3.0  kcal  mol.  falls  in  the  2-5  kcal  mol  range  predicted  [9]  and  confirmed  by 
theoretical  studies  [17-20].  This  barrier  is  low  enough  to  suggest  rapid  interconversion  between 
the  two  preferred  sugar  conformations  but  high  enough  to  imply  that  pseudorotation  is  hindered. 
The  barrier  observed  at  the  West  is  significantly  higher  and  suggests  that  conformational 
interchanges  through  the  unfavorable  Ol'-exo  conformation  are  prohibited. 


ENERGIES.  \n  examination  of  the  energy  component  curves  can  explain  the  observed 
trends  in  the  relative  energy  changes  for  the  8  different  model  sets  studied  (see  Figure  7). 
Clearly,  the  overall  shape  of  the  total  energy  curve  is  dominated  by  that  of  the  non-bonded 
component.  However,  if  only  the  nonbonded  interactions  are  considered,  the  energetically  equal 
C2'  — endo  and  C3'  — endo  conformations  are  separated  by  a  very  low  East  barrier  which 
consequently  suggests  that  pseudorotation  is  free.  Only  the  bond  strain  and  threefold  torsional 
potential  can    effectively  raise  the  East  and  West  barrier  heights. 

Both  the  bond  length  and  bond  angle  energies  exhibit  sinusoidal  dependencies  as  expected 
from  the  empirical  functions  used  to  describe  the  variations  K^i  these  geometric  quantities  with  P 
[66].  .As  noted  previously  [19].  the  threefold  torsional  potential  constitutes  the  largest  contribution 
in  magnitude  to  the  total  energy  along  the  pseudorotation  pathway.  The  sinusoidal  torsional 
potential  exhibits  2  local  minima  at  O'  and  180'  .  corresponding  to  the  symmetrical  Twist  forms  of 
C2"-endo  -  C3'  — exo  and  C3'  — endo  -  C2"  — exo  (see  Fisure  3). 
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The  gauche  potential  contributes  0.*^5  kcal  mol  at  P  =  0*^'.  increases  slowly  until  2"'''.  then 
falls  to  0.1  kcalmol  at  164'  ,  and  rises  back  to  0.95  kcal  mol  to  complete  one  pseudorotation  cycle. 
Thus,  the  energy  difference  of  the  gauche  contribution  between  the  North  and  South  conformers  is 
essentially  the  same  as  that  observed  by  minimization  -  about  0.85  kcal  mt)l. 

The  effects  on  the  energy  differences  between  the  two  local  minima  can  now  be  understood. 
From  Table  IV  it  is  evident  that  the  energy  difference  between  the  C2"-endo  and  C3'-endo 
structures  is  affected  by  several  changes  in  the  energy  function  form  and  as>ociated  parameters. 
Introduction  of  a  gauche  potential,  a  weighted  threefold  torsional  potential,  or  a  lowering  oi  the 
bond  stretching  constant  can  effectively  lower  the  energ_\  of  the  C2'  — endo  conformation  in 
relation  to  C3"-endo. 

Since  the  bond  energy  is  sinusoidal  and  non-negative  with  a  minimum  of  0.0  kcalmol  at 
P  =  17*  .  a  lowering  of  the  stretching  parameter  S\  hardly  affects  the  C3"-endo  energy.  The 
function's  sinusoidal  shape  is  dampened,  thereby  lowering  the  energy  for  C2'-endo.  This 
explains  the  energy  shift  for  sets  6.  7  and  8  of  the  E  vs.  P  study. 

By  modeling  all  dihedral  angles  and  appropriately  adjusting  the  torsional  parameters  \'}.  the 
torsional  curve  is  essentially  translated  downward  by  about  1.2  kcal  mol.  Fluctuations  in  this  shift 
are  simply  due  to  the  difference  in  the  overall  dihedral  angle  contribution  for  various  P  values. 
Since  the  North  minimum  is  lowered  by  1.08  kcal  mol  and  the  South  minimum  by  1.27  kcal  mol. 
the  net  affect  of  this  torsional  potential  modification  is  to  slightly  lower  the  energy  at  the  South. 
The  observed  translation  of  the  E-i.fQi[)  curve  indicates  that  equal  treatment  of  all  heavy  atom 
rotations  can  be  achieved  correctly' by  this  weighing  procedure. 

Finally,  we  can  explain  the  observation  that  changing  the  bending  constant  .S'2  has  no  effect 
on  the  energy  difference.  By  increasing  .V2  from  15.0  to  hO.O  kcalmol  for  N'-C-C.  N'-C-O. 
and  all  bond  angles  involving  hydrogens,  the  biMid  angle  bending  cnergv  is  increased  uniformlv  by 
about  0.4  kcalmol.  Thus,  the  energy  difference  between  the  two  pret'erred  ci)nforinations  is  left 
unchanged. 


-  3Q  - 
4.4     Comparison  of  the  Two  Approaches 

The  energy  differences  between  the  C3'  — endo  and  C2"-endo  conformations  are  not  identical 
for  the  energy  minimization  and  E  vs.  P  studies,  but  they  exhibit  the  same  relative  energy 
difference  trends  i  see  Table  IV  i.  For  both  studies,  the  gauche  potential  can  be  used  to  shift  the 
global  minimum  from  the  C3"  — endo  to  the  C2"  — endo  structure  by  0.85  kcal  mol.  .A  threefold 
torsional  potential  for  all  heavy  atom  dihedral  angles  can  also  shift  the  global  minimum  to  the 
South,  by  about  0.2  kcal  mol.  However,  since  the  energy  differences  are  larger  for  the  E  vs.  P 
curves  than  for  the  minimization,  the  energy  shift  induced  b\  the  weighted  threefold  term  makes 
the  two  minima  energetically  equal  in  the  E  \s.  P  curves. 

In  contrast  to  the  minimization  study,  changing  the  bond  stretching  parameter  ."^1  affects  the 
energy  difference  between  the  two  local  minima  in  the  E  vs.  P  curves.  .As  ue  see  in  ^ets  -i-  and  6 
in  Table  IV.  lowering  .S' 1  from  100.0  to  25.0  kcal  mol  .^  ^  hardly  changes  the  energy  of  the 
minimization-generated  structures  but  it  effectively  raises  the  energy  difference  by  about  0.23 
kcal  mol  for  the  E  vs.  P  results.  Similarly,  the  relative  energy  changes  not  observed  in  the 
minimization  <ets  5  and  S  are  contrasted  by  an  increase  of  about  0.2  kcal  mol  in  favor  of  the  South 
minimum. 

Thus,  these  observations  demonstrate  that  in  general  the  Energy  vs.  P  approach  is  more 
susceptible  to  fine-tuning  and  manipulation  for  reproducing  the  desired  quantitative  energetic  and 
structural  data.  .Any  geometric  differences  between  the  structures  along  the  pseudorotation  cycle 
can  be  exploited  by  appropriate  choices  of  energy  components  and  parameters.  This  would  also 
be  the  case  v\ith  minimized  E  vs.  P  curves,  since  geometry  is  generated  by  constraints.  In  energy 
minimization  with  the  full  set  of  degrees  of  freedom,  the  molecule  can  effectiveK  perform 
"internal  adjustments"  to  accomodate  the  change^  in  energy  parameters. 

In  summary,  results  from  minimization  and  E  vs.  P  profiles  indicate  that  the  use  of  a  gauche 
potential  and  a  weighted  threefold  torsional  potential  are  appropriate  for  reproducing  an 
energetically  more  favorable  state  for  deoxynbose  at  the  South  region  of  the  pseudorotation  cycle 
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and  tor  modeling  systematically  all  heavy  dihedral  rotations.  Thus,  parameter  sets  5  or  8  produce 
deoxyribose  conformations  that  are  consistent  v\ith  structural,  energetic  and  geometric  data 
observed  in  crvstal  structures. 


4.    CONCLUSIONS 

Previous  discrepancies  witin  experiment  for  the  furanose  ring  conformation  have  prompted 
our  present  investigation.  These  discrepancies  have  revealed  some  underl\ing  uncertainties  ot 
potential  energy  calculations.  Thus,  we  have  chosen  to  examine  the  fundamental  issues  ot  a 
computational  approach  -  degrees  of  freedom,  potential  energy  parameterization  and 
minimization  for  the  deo\yribo>;e  model.  The  detailed  examination  ot  >tructural.  energetic  and 
geometric  results  for  various  energy  parameter  sets  b>  minimization  and  E  \>.  P  cur\ev  has 
allowed  us  to  evaluate  the  pseudorotation  approximation  and  the  influence  of  various  parameter 
and  modeling  strategies. 

For  any  computational  approach,  a  reduction  of  the  full  set  of  3.V-6  degrees  of  rreedom 
associated  with  a  molecule  of  .V  atoms  involves  a  trade-off  between  accuracy  and  reliability  of 
potential  energy  studies  on  one  hand  and  computational  simplicity  on  the  other.  For  nucleic  acid 
sugars,  a  reduction  accomplished  by  the  pseudorotation  approximation  has  apparently  provided  a 
good  qualitative  description  of  various  puckering  modes.  However,  this  simplified  analytic 
formulation  has  undeniable  drawbacks. 

First,  fluctuations  m  puckering  amplitudes  have  been  theoreticall>  and  experimentally 
observed  [6.7.19.34].  Second,  the  analytic  development  of  pseudorotation  is  based  on  statistical 
analvsis  of  susav  fragments  in  nucleic  acids;  various  averaging  techniques  are  used  to  fit  structural 
data  (bond  lengths,  bond  angles  and  dihedral  angles)  of  many  compounds  to  analytic  expressions 
v\ith  the  same  coefficients.  Third,  extension  of  the  pseudorotation  approximation  to  large 
molecular  systems  is  difficult  and  compoundedly  approximate.  The  inclusion  o(  P  and  possibly 
~ma.\  ^^  independent  variables  of  the  energy  function  requires  derivatives  to  be  computed  u  ith 
respect  to  P  (and  ~mM  'f  't  '^  variable)  in  each  minimization  step.  If  the  conformational  energy  is 
described  in  cartesian  coordinate  space,  this  implies  a  complicated  use  of  the  chain-rule. 
.Alternatively,  if  an  internal  coordinate  system  of  dihedral  angles  is  chosen,  each  value  ot  P  and 
Tmjv  must  be  used  as  conformational  variables  in  some  coordinate-generation  procedure  tor  all  the 


atoms  in  the  molecule.  For  all  these  reasons,  the  sacrifice  of  3  degrees  of  freedom  (for  the  5 
dihedral  angles i  hardly  seems  justified.  If  analysis  of  >tructural  data  for  various  sugar  fragments 
is  desired,  the  e.xact  Fourier  series  ill)  can  provide  a  simple  but  more  reliable  .scheme. 


With  regard  to  the  potential  energy  construction  and  energy  minimization,  the  overall 
structural,  energetic,  and  geometric  agreement  of  our  minimized  structures  with  experimental  data 
suggests  that  the  potential  energy  function  is  well  constructed  and  that  the  Newton  methods  are 
powerful  minimization  techniques.  By  a  modeling  strategy  that  involves  a  simple  moditication  to 
the  usual  harmonic  bond  length  and  bond  angle  potentials,  we  construct  expressions  that  are 
simpler  and  hence  easier  to  evaluate.  This  in  turn  reduces  the  potential  for  round-off  error.  To 
eliminate  the  arbitrariness  in  modeling  dihedral  angles,  v^e  ha\e  implemented  equal  treatment  of 
all  heavy  rotations  in  the  molecule  and  have  demonstrated,  using  the  E  \s.  P  cur\es.  that  such  a 
treatment  is  indeed  valid.  In  combination  with  polynomial  representations  of  the  bond  angle 
cosine  and  dihedral  angle  cosine  expressions  in  terms  of  the  coordinate  \ariables  and  with  group 
treatment  of  bond  stretching,  angle  bending,  and  dihedral  angle  parameters,  we  avoid  generation 
of  unrealistic  geometry  and  make  feasible  an  application  of  an  efficient  minimization  strategy 
which  uses  analytic  second-derivatives.  We  are  currently  examining  these  modeling  and 
minimization  issues  for  larger  DNA  seaments. 
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FIGURE  LEGENDS 


Figure  1.  A  schematic  representation  of  the  2'- deoxynbose  model  investigated  in  this  study ^  The  base 
of  a  nucleoside  is  replaced  by  an  amino  group,  and  the  CH^OH  group  attached  to  the  C4'  atom  of  a 
nucleoside  is  replaced  by  a  methyl  group.  The  sugar  atoms  are  designated  by  primes  as  in  the  conventional 
nucleic  acid  notation  to  distinguish  sugar  atoms  from  base  atoms.  9,^  -  9_  are  the  5  endocyclic  bond  angles  at 
atoms  or  through  C4',  in  clockwise  order,  and  t,i  -  t_  are  the  5  endocyclic  dihedral  angles  defining  the 
rotations  about  bonds  Ol'-Cl',  Cl'-CZ',  C2'-C3'.  C3"-C4'  and  C4'-Or  .  respectively.  .A  dihedral  angle 
is  defined  0"  if  all  4  atoms  involved  are  coplanar.  The  sign  of  the  dihedral  angle  is  determined  by  the 
procedure  described  in  Section  3.1. 

Figure  2.  The  two  common  sugar  puckering  modes  of  C3'-endo  and  C2"-endo  as  generated  by 
minimization.  Two  different  views  are  given  for  each  conformation.  In  classifying  puckered  conformations 
for  five-membered  rings,  two  basic  forms  are  discussed:  an  envelope  (E)  in  which  tour  atoms  are  in  a  plane 
and  the  fifth  atom  out.  or  a  twist  (T).  where  two  adjacent  atoms  are  displaced  on  opposite  sides  from  the 
plane  defined  by  the  other  three  ring  atoms.  Sugar  puckering  modes  are  defined  according  to  which  atom  or 
atoms  are  displaced  from  these  three  or  four  atom  planes:  atoms  displaced  on  the  same  side  as  e.xocyclic 
carbon  C5'  are  called  endo ,  and  those  on  the  opposite  side  of  Cf  are  termed  exo. 

Figure  3.    The  pseudorotation  cycle  of  the  furanose  ring.    Ten  symmetrical  Twists  (T)  are  defined  for 

phase    angle    values    .of    P  =      0",36'\72'' 324".    and    ten    Envelopes    (E)    are    defined    for    P   = 

18'',54'\90'^ 342'^.     Twelve    illustrative    puckering    forms    are    drawn    -    all    10    Envelopes    and    the    2 

symmetrical  Twists  of  C3'- endo- C2'-e.xo  at     P  =  0"  and  C2'- endo- C3'- exo  at     P  =   180". 

Figure  4.  Definition  of  a  dihedral  angle.  The  dihedral  angle  t,,;,  is  defined  as  the  angle  between  the 
normal  to  the  plane  defined  by  atoms  i-j-k  and  the  normal  to  the  plane  defined  by  atoms  y- A-/. 

Figure  5.  A  schematic  representation  of  positioning  the  exocyclic  substituents  of  ring  atoms.  Atom  /. 
attached  to  ring  atom  j,  is  positioned  in  a  local  bisecting  coordinate  system  (e,  .e,  .e-J  through  the  3  ring 
atoms  ;  -  7  -  i.  as  follows.  Let  a  =  .\,  -  x.  ,b  =  x,  -  x  and  r  be  the  distance  between  x.  and  x,.  Then 
we  have: 

„    -  (a  -   b)  (a  X   b) 

We  define  (j),  as  the  angle  between  p  and  the  e,  axis,  where  p  is  the  projection  vector  of  x  with  the  e,  ,  e, 
plane;  <1),  is  defined  as  the  angle  between  x  and  p  .  Once  <b ,  and  <J),  are  specified,  the  coordinates  of  x  are 
given  by: 

X/  =  x^  -   (r    cos4),  sin<|>|)  e,   -   (r  .  cos<j),  coscj),)  e,  -  (r,  sindO  ^\  ■ 

To  position  the-  exocyclic  atoms  symmetrically  from  the  ring  plane  so  that  a  tetrahedral  arrangement  is 
approximately  occupied,  we  let  (b,  -  0  and  <}),  ^  r  54.7356.  For  C5',  N  .  H3'  and  1H2'  *,  -  -  54.7356 
and  for  H4',  HI',  03'  and  2H2'  d),  -   -  54.7356(see  Figure  1). 
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Figure  6.  Illustration  of  the  minima  obtained  from  starting  points  along  the  pseudorotation  cycle. 
Startmg  pomts  at  phase  angle  values  at  18"  intervals  generate  the  nearest  local  minimum  on  the 
pseudorotation  path. 

Figure  7.  Results  for  Energy  vs.  P  curves.  The  total  potential  energy  E  and  its  components  are 
illustrated  as  a  function  of  the  phase  angle  of  pseudorotation  P  for  the  model  deoxyribose  studied.  Results 
are  given  from  6  variations  of  the  potential  energy  function  corresponding  to  sets  2,3,4,5,6  and  8  of  Table  II. 


TABLE    LEGENDS 


Table  I.  Dihedral  angles  and  associated  parameters  investigated  in  this  study.  In  column  fb)  of  the  V3 
parameter,  all  heavy  atom  rotations  in  the  model  are  included,  with  V3  appropriately  scaled  from  values  in 
column  (a),  where  only  endocyclic  rotations  about  C-C  and  C-O  bonds  are  included.  All  angles  are  given 
in  degrees  and  all  energy  parameters  in  kcal/mol. 

Table  11.  A  description  of  energy  functions  used  in  this  study.  The  coefficients  for  the  non-bonded 
interactions  remained  constant  throughout  the  calculations,  and  variations  in  bond  length,  bond  angle  and 
torsional  terms  comprise  the  various  sets.  The  bond  length  force  constant  51  is  given  in  kcal/mol  A^  and  the 
bending  and  torsional  parameters  in  kcal/mol.  The  angle  bending  force  constant  52  has  a  value  of  60 
kcal/mol  for  all  C-O-C,  C^C-O  and  C-C-C  sequences  throughout  the  sets.  Only  the  bending  constants 
for  N-C-O  and  N-C-C  angles,  and  all  sequences  involving  hydrogens  (denoted  by  H**)  were  varied. 
The  references  in  the  torsional  columns  correspond  to  the  appropriate  columns  of  Table  I. 

Table  III.  A  summary  of  geometric  and  energetic  results  from  the  minimization  study.  For  all  8 
energy  sets  described  in  Table  II.  various  quantities  associated  with  the  generated  North  (N)  and  South  (S) 
minima  are  given:  the  estimated  phase  angle  of  pseudorotation  P  and  the  puckering  amplitude  t^^.^ 
(obtained  from  the  procedure  oulined  in  Section  4.1).  the  endocyclic  bond  angles,  endocyclic  dihedral  angles. 
2  dihedral  angles  associated  with  rotations  about  the  C3'-C4'  bond,  and  all  energy  components.  For 
comparison,  we  report  bond  angles  obtained  experimentally  for  the  C3'-endo  and  C2'-endo  conformations 
[7.17],  and  endocyclic  dihedral  angles  predicted  by  the  pseudorotation  equation  (2)  with  7,^^^  =  38''  .  All 
angles  are  given  in  degrees  and  all  energies  in  kcal/mol. 

Table  IV.  A  comparison  between  the  minimization  and  E  vs.  P  studies.  For  all  8  energy  sets 
examined,  the  values  of  P  and  E  are  given  for  the  N  and  S  minima  as  well  as  the  energy  difference. 
£y  -  E^,  for  each  set.    All  angles  are  given  in  degrees  and  all  energies  in  kcal/mol. 
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Dihedral  Angles  and  .Associated  Parameters 


;1  Twofold  Torsion  Constants  (V2) 

Threefold  lorsion  Constants  (\'3) 

Dihedral  .Angle 

'i 

a                                 b                           1 

(-,)     C4-or-cr-c:2' 

ij 

1.8                             0.9 

C4-Or-Cl-N 

1 

0.9 

1 

(t,)      oi-cr-C2'-C3- 

;                            0- 

2.8                              1,4 

N-Cr-C2-C3' 

1 

1.4 

(T,)       Cl-C2'-C3--C4- 

1! 

2.8                                1.4 

Cr-C2-C3-03- 

0.2 

II 

l! 

1.4 

(7,)       C2'-C3--C4--Or 

i|                                    0.2 

2.8                               0.7 

C2-C3'-C4-C5' 

1 
i 

0.7 

03--C3-C4-Or 

:,                             1.0 

0.7 

03--C3--C4--C5- 

1 

0  2 

0.^                        1 

1 

(tj     c3--C4'-or-cr 

1! 
il 

1.8                             0.9 

c5-C4'-oi  -cr 

,i 
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"                                                                              ! 
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Descripdoa  of  Energy  Sets 


SET 

Paramcteri 
Modified 

SI 
(Bond  Length) 

S2 

(Bond  Angle) 

V2 

(sec  table  V) 

V3 

(see  table  V) 

V3.a 

1 

100  for  ail  bonds 

15  for  NCO  and  NCC 
200  for  H" 

2 

S2 

. 

15  for  NCO,  NCC 

and  H" 

3 

V2 

•• 

- 

V2,a 

- 

4 

V3 

M 

m 

V3.b 

5 

V2.V3 

- 

- 

V2,a 

- 

6 

S1.V3 

25  for  all  bonds 

1* 

- 

i. 

7 

S1.S2.V3 

- 

60  for  all  angles 

- 

" 

8 

S1,S2,V2,V3 

50  for  all  bonds 

- 

V2,a 

" 

Table    III 

Mininiizarion  Results  -    Geometric  and  Energetic  Details 


SET 

Experimental 

1 

2 

3 
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Minimum 

C3 '-endo 

C2'-endo 

N 

S 

N 
'     17 

S 

N 

S 

N 

24 

S 

? 

18 

162 

18 

168 

162 

9 

162 

1   158 

•q 

110 

110 

110 

111 

110 

110 

110 

110 

109 

109 

9, 

107 

106 

105 

104 

106 

105 

106 

105 

107 

105 

9, 

102 

101 

104 

103 

102 

102 

102 

102 

102 

102 

e, 

102 

103 

101 

102 

101 

102 

101 

102 

101 

103 

«4 

104 

106 

105 

107 

105 

107 

106 

107 

104 

107 

Tn 

0 

-22 

0 

-20 

-4 

-23 

7 

-23 

-4 

-25 

T| 

-22 

36 

-24 

35 

-28 

37 

-29 

37 

-22 

37 

T-, 

36 

-36 

37 

-37 

39 

-36 

39 

-36 

37 

-35 

''y 

-36 

22 

-37 

26 

-38 

23 

-36 

23 

-40 

21 

T« 

22 

0 

24 

-4 

22 

-1 

19 

0 

28 

2 

mi» 

38 

38 

39 

38 

41 

155 

39 

40 

39 

41 

3S 

01'-C4.C3-03' 

1 

157 

95 

96 

153 

96 

159 

96        1 

c«-C4'-a'.03' 

83 

146 

85 

143 

88 

144 

82 

146 

Esn 

-0.34 

0.31 

-1.11 

-0.87 

-1.08 

-0.86 

-1.27 

-0.74 

^ao\D 

0.09 

0.11 

0.04 

0.04 

0.04 

0.04 

0.06 

0.07 

^BA.\C 

4.20 

4.03 

4.17 

3.77 

4.11 

3.81 

4.16 

3.52 

^IFQUi 

0 

0 

0 

0 

1.30 

0.46 

0 

0 

^iFOLD 

8.13 

8.27 

7.77 

8.32 

11.27 

7.82 
12.20 

8.27 
11.73 

6.86 

9.81 

7.21 
10  06 

^  TOTAL 

12.08 

12.73 

10.87 

Pieudoroimon 
E>eviaaons 


IM 
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Tabic  III  (coDtinued) 

MinimizitioD  Resulti  •   Geometric  and  Energetic  Details 


SET 

EipenmentaJ 

5 

6 

7 

8 

Mmunum 

C3'-endo 

C2"-endo 

N 
20 

S 

N 

S 

N 
20 

S 

N 

S 

P 

18 

162 

159 

25 
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18 
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«. 

no 
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101 
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101 
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9. 
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e. 
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Tn 

0 

-22 

-1 

-25 

-4 

-26 

-1 

-26 

-1 
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^: 
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36 
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37 

-21 

38 

-23 

39 

-25 

38 

T, 

36 

-36 

37 

-35 

37 

-36 

37 

-36 

37 

-36 

T] 

-36 

22 

-39 

22 

-41 

21 

-39 

22 

-37 

23 

T4 

22 

-  ..  - 

0 

25 

2 

28 

3 

25 

2 

23 

1 

mi» 

38 

38 

40 

39 

41 
159 

35 

40 

40 

39 

40 

or-C4-.a'-03' 

157 

96 

96 

158 

96 

157 

95 

a-C4'-a-o?- 

83 

146 

81 

146 

81 

146 

83 

146 

f'.va 

-1.25 

-0.73 

-1.41 

-0.92 

-0.30 

0.18 

-0.16 

-0.34 

^BOSD 

0.05 

0.07 

0.22 

0.27 

0.23 

0  29 

0.10 

0.15 

^B\.\C 

4.07 

3.54 

4.10 

3.44 

4.10 

3.73 

4.04 

3.77 

^  IF  OLD 

1.37 

0.48 

0 

0 

0 

0 

1.03  1 

0.14 

^irOLD 

6.96 

7.18 

6.73 

7.06 

6.94 

7.00 
11.20 

7.12 
12.14 

7.09 
11.49 

^  TOTAL 

11.19 

10.54 

9.63 

9.85 

10.97 

Pjeudorotiaon 
Deviaaoni 

1.20 

2.58            1-21 
1 

59.99 

1.20 

2  29 

4  70 

2.67 

Table     IV 


Summary  of  Minimization  and  Energy  vs.  P  Results 
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Parameters 
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Lo 
£  ^ 
N 
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^s.  ? 
S 

of  Minima 
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N          S 

£  vs.  P 

MinimLzaaon 

1 

- 

15 

177 

18       168 

-0.21 

-0.65 

2 

52 

14 

176 

17       162 

-0.15 

-0.40 

3 

V2 

13 

176 

9       162 

*0.72 

-0.53 

4 

V3 

15 

174 

24 

158 

-0.01 

-0.25 

5 

V2,V3 

14 

174 

20 

159 

♦0.89 

-0.65 

6 

S1,V3 

14 

166 

25 

162 

-0.34 

-0.22 

7 

S1,S2,V3 

14 

166 

20 

158 

-0.32 

-0.23 

8 

S1.S2,V2,V3 

14 

169 

18       160 

1.06 

-0.65 
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APPENDIX  A 
GENERAL  CONCEPTS  OF  NONLINEAR  OPTIMIZATION  METHODS 


In  this  section,  we  attempt  to  present  only  the  key  ideas  in  minimization  methods  and  categories.  For 
general  references,  the  reader  can  consult  [39-44|.  In  these  books,  especially  the  more  recent  ones,  many 
references  for  specific  algorithms  that  are  mentioned  below  can  be  found. 

(a)     Mathematical  Preliminaries 


In  the  following  discussion  we  generally  denote  scalars  by  lower-case  Greek  letter  such  as  a  ,  3  .  \  : 
column  vectors  by  bold  faced  lower-case  Roman  letters  as  .x  ,  p  ,  and  matrices  by  capital  Roman  letters.  .A 
superscript  T  denotes  a  vector  transpose,  so  that  x'  represents  a  row  vector  and  x  y  an  inner  product.  .A 
superscript  T  also  denotes  a  matri.x  transpose. 

Unless  stated  otherwise,  all  vectors  belong  to  R"   ,  the  n-dimensional  vector  space.  The  standard  basis 

vectors   in   R"   are   the  n   vectors   {  e,  ,  e, e,  }   where  e    is  the  vector  with  the  entry    1   in   the  ;-th 

component  and   0   in   all   others.    The   vector  norm    that  we   associate   with   a  vector  in   R"    is   the   standard 


Euclidian    norm,  which  we  denote  bv 


and  define  as    x 


Ik  I- 


We  are  interested  in  solving  the  optimization  problem  without  constraints 


minimize      E(\)  .  x  €  R" 


(1) 


locally,  where  £  is  a  real-valued  function  of  n  variables:  E(\)  -  E(-r,  ,  .t,  ,  ...  .r.,).  We  will  assume  in  the 
following  discussion  that  E  is  sufficiently  smooth,  that  is.  E  possesses  continuous  derivatives  up  to  order  m 
where  m  >  3.     The  gradieni   of  £   at  x     is   defined     to   be   the   first-derivative  vector  V£(x)     or    g(x)  of  n 

components       .g,(-x)  =   ^^ — -,  and   the  Hessian  at      x   is  defined   to   be   the  n\n   mairi.x  of  second  partial 


dx. 


derivatives    denoted    as     H(x)    with    components       W   (x) 
Hessian  matri.x  is  symmetric ,  i.e.  H   (\)   =  //,,(x). 


(5.r    (?.r. 


Since 


a-£(x)   „   a-£(x) 


a.r,  a.r. 


dx.  a.t. 


the 


A  point  X  6  R"  is  said  to  be  a  strict  local  minimum  of  £  if  there  exists  an  ti  >  0  such  that 
£(x")  <  £(x)  for  all  xi*:  x'  within  a  distance  ti  of  x  .  i.e.  ||x  -  x"  ||  <  ti.  The  point  x'  is  a  global  minimum  of 
£  if  £(x')  <  £{x)  for  all  xT^x",  where  x  h   R  '. 


There  are  several  general  characteristics  of  a  matrix  that  are  particularly  useful  for  analysis  of 
minimization  algorithms. 

Density  of  a  matrix  is  a  measurement  given  by  the  ratio  of  the  nonzero  to  zero  matrix  components.  .A 
matrix  is  said  to  be  dense  when  this  ratio  is  large  and  sparse  when  it  is  small. 

A  symmetric  matrix  A  is  said  to  be  positive -dejiniie  if  the  quadratic  form  u'^.Au  >  0  for  all  nonzero 
vectors  u.  Similarly,  A  is  negative-definite  if  u'^Au  <  0  for  all  nonzero  vector  u.  A  is  indefinite  if  u^.Au  is 
positive  for  some  u  and  negative  for  others. 
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The  Hessian  matrix  is  a  generalization  to  R''  of  the  concept  of  curvature  of  a  function.  The  positive- 
definiteness  of  the  Hessian  is  a  generalized  notion  of  positive  curvature. 

We  now  wish  to  characterize  the  sufficient  conditions  that  must  hold  at  a  solution  point  of  (1).  These 
conditions  are  simply  extensions  to  R'  of  the  well-known  first  and  second-derivative  conditions  for  a 
univariate  function. 


Theorem:    Let  E{\)  ,  x  €  R"  be  a  smooth  function  with  continuous  first,  second  and  third-derivatives 
defined  for  all  x.  Suppose  the  following  conditions  hold  for  the  gradient  and  Hessian  of  E  at  \'  ■;  R": 

(0      g(x")  =  0 

(u)     H{\')    IS  positive-definite. 
Then  x    is  a  strict  local  minimum    of  E. 

Proof:    By  Taylor's  Theorem  with  remainder,  we  can  expand  E(\)  along  an  arbitrary  direction  p  f  R": 
£{x"-pl  =   £(x-)  -  g(x   )'p  -   :J-p"W(x-|p  -  0(|!pjp)  .  (2) 

■'?',">■■ 

Since  condition  (/)  holds,  the  gradient  term  vanishes,  and  condition  (//)  implies  that  there  exists  some  a  >  0 
such  that  p^/y(x")p  >  a||p||-.    Then  the  Taylor  series  can  be  written  as: 

£(x--p)  -  Eix)  >   f  llpiP  -   0(||pip)  . 

For  arbitrary  small  ||p||  .  the  first  term  on  the  right  dominates  the  second,  and  consequently  the  right  hand 
side  is  positive.  By  definition,  x'  is  a  strict  local  minimum  of  E . 

In  analogy  to  the  one-dimensional  case,  a  point  \     is  a  siatwnarx  point  of  E  if  g(x  i   =   0  but  H {\  )  is 
not  positive-definite. 


(b)     Basic  Descent  .Methods 

The  fundamental  structure  for  iterative  techniques  for  solving  unconstrained  minimization  problems  is 
simple.  A  starting  point  is  chosen:  a  direction  of  movement  is  chosen  according  to  some  rule,  and  a 
line  search  is  performed  to  approximately  minimize  the  objective  function  along  that  direction.  The  process 
is  repeated  at  the  new  point,  and  the  algorithm  continues  until  a  local  minimum  is  found.  Schematically,  a 
model  descent  method,  can  be  written  as  follows: 

1.  Supply  an  initial  guess  x,, 

2.  For  i  =  0,  I,  2,  .  .  .  until  convergence 

*  determine  a  search  direction  p,  V-^/ 

*  determine  an  appropriate  step  length  \,. 

*  set   x,.,,   =   xj  -  \.  p,. 


Algorithms  differ  by  their  choice  of  search  directions  and  consequently  have  varying  computational  costs  and 
convergence  properties. 

VIost   algorithms   are    analyzed   by    investigating   their   application   to   a   quadratic  function.    A   general 
quadratic  function  m  multivariate  space  can  be  written  as 

E{\)  =   -  \'A\  -  \yx  (4) 


where  A  is  a  positive  definite  n^n  matrix,  b  a  constant  vector  in  R".  and  x  a  general  vector  in  R".  The 
gradient  of  this  function  is  the  vector  \\  -  b.  and  the  Hessian  is  the  matrix  A.  According  to  the  Theorem  of 
the  previous  section,  a  local  minimum  of  E(\)  can  be  found  by  setting  the  gradient  equal  to  zero.  Such  a 
point  X    will  then  satisfy  the  linear  system: 

Ax"   =   -  b  (5) 

For  positive-definite  matrices  .4.  this  solution  is  unique,  and  consequently  ,\  is  the  global  minimum  of  £. 
Since  general  functions  are  well  approximated  by  a  quadratic  function  in  the  neighborhood  of  their  local 
minima,  the  convergence  properties  obtained  for  quadratic  functions  can  be  extended  to  general  functions 
locally. 

The  convergence  properties  of  an  algorithm  are  described  by  two  analytic  quantities:  convergence  order 
and  convergence  ratio.  Let  the  sequence  x,-,  ,  x,  .  x.  ,  .  .  .  converge  to  x".  i.e.  lim  ..,.  ||x;.  -  x'|i  =0.  The 
sequence  is  said  to  converge  to  x"  with  order  p  if  p  is  the  largest  non-negative  number  for  which  the  limit  3 
exists  and  is  finite;  where 

0  s  3   =  lim      -,^—^ .  (6) 

■""iix, -x-ir 

When  p  =  1  and  3<  1  the  sequence  is  said  to  converge  linearly  (e.g.  x.=2~'./i  =  l),whenp  =  land  3=0 
the  sequence  converges  superlinearly  (e.g.  x.  =i~'  .  n  =  l).  and  when  p  =2  the  convergence  is  quadratic 
(e.g.    X;=2~-    ,  :i==  1).    The  constant     3  '^  the  associated  convergence  ratio. 


(c)     Non-derivative  Methods 


Vtinimization  methods  that  incorporate  only  function  values  generally  involve  some  systematic  method 
to  search  the  conformational  space.  In  Coordinate  Descent  methods,  the  search  directions  are  the  standard 
basis  vectors.  A  sweep  through  these  n  search  vectors  produces  a  sequential  modification  of  one  function 
variable  at  a  time.  By  repeatedly  sweeping  the  n-dimensional  space,  a  local  minimum  might  ultimately  be 
found. 

A  well-known  variant  of  the  basic  coordinate  descent  scheme  is  Powell's  method.  The  standard  basis 
directions  are  modified  as  the  algorithm  progresses  so  that  when  the  procedure  is  applied  to  a  quadratic 
function,  it  generates  n  mutually  conjugate  directions  after  n  sweeps.  A  set  of  mutually  conjugate  directions 
possesses  the  important  property  that  a  search  in  each  of  the  directions  only  once  suffices  to  find  the 
minimum  solution.   Powell's  method   thus    guarantees  that  in  exact  arithmetic  (i.e.   in  absence  of  round-off 


error)  the  minimum  of  a  quadratic  is  found  after  n  sweeps. 

Minimization  methods  in  this  class  are  generally  easy  to  implement  and  do  not  require  any  derivative 
computation,  but  their  realized  convergence  properties  are  rather  poor.  They  may  work  well  in  special 
cases  when  the  function  is  quite  random  in  character  or  the  variables  are  essentially  uncorrolated.  In 
general.  The  computational  cost,  dominated  by  the  number  of  function  evaluations,  can  be  excessively  large 
for  functions  of  many  variables  and  can  far  outweigh  the  benefit  of  avoiding  derivative  calculations. 

If  obtaining  the  analytic  derivatives  is  out  of  the  question,  viable  alternatives  still  remain.  The  gradient 
can  be  approximated  by  finite -differences  and  used  as  such  in  a  gradient  or  quasi-Newton  method.  These 
alternatives  will  generally  provide  significant  improvement  in  the  computational  cost,  and  also  lead  to  local 
minima  that  might  otherwise  remain  undetected.  Despite  all  the  drawbacks  of  non-derivative  methods,  their 
ease  of  application  has  made  them  the  choice  for  some  potential  energy  applications. 


Id)      Gradient  Methods 

Methods  that  use  analytic  values  for  derivatives  clearly  possess  more  accurate  information  about  the 
function.  Gradient  methods  can  use  the  slope  of  a  function  as  directions  of  movement  toward  cxtrcmum 
points;  second  derivative  methods  can  use  curvature  information  directly  from  the  Hessian  in  order  to  find 
the  region  where  the  function  is  convex.  Two  common  gradient  methods  are  Steepest  Descent  (SD)  and 
Conjugate  Gradient  (CG). 

Steepest  Descent.  Steepest  Descent  is  one  of  the  oldest  and  sin)plest  methods.  It  is  actually  more 
important  theoretically  than  practically  as  a  reference  by  which  to  test  other  methods.  At  each  iteration  of 
Steepest  Descent,  the  search  direction  is  taken  as  the  negative  gradient  of  the  objective  function  at  X; ,  -  g;.. 
To  understand  this  choice,  consider  the  linear  approximation  to  E  at  X;  based  on  a  truncated  Taylor  series 
along  a  search  direction  p  : 


E(x_ 


£(x,-p)  =  E(x,)  -  g.  ^ 


(7) 


It  is  reasonable  to  choose  a  search  direction  p:   that  will  be  a  descent  direction  in  the  sense  that  the  new  x 
iterate  will  cause  a  reduction  in  the  function  value.    This  is  equivalent  to  the  condition  that 


g.  '  P  <  0 


(8) 


The  simplest  way  to  guarantee  this  is  to  choose 


Pi 


(9) 


This  choice  minimizes  the  inner  product    -  g,^  p  for  vectors  of  unit  length,  and  hence  the  name  steepest- 
descent. 

SD  is  simple  to  implement,  requires  little  storage  (order  of  n).  but  may  be  very  slow,  especially  near  a 
solution.    The  convergence  rate  of  SD  when  applied  to  a  quadratic  function  is  only  Imeur  with  a  convergence 

where  r,  the  condition  number,  is  the  ratio  of  largest  to  smallest  eigenvalues  of 


ratio  no  greater  than 
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A  :  r  =  \„,  .J.4)  /  \^,,„(A).  When  r  is  large,  the  convergence  ratio,  which  measures  the  reduction  of  error  at 
every  step,  can  be  arbitrarily  close  to  1.  In  other  words,  the  convergence  rate  of  SD  is  slowed  as  the 
contours  of  the  obejective  function  become  more  eccentric. 

Conjugate  Gradient.  The  Conjugate  Gradient  method  was  originally  designed  to  minimize  quadratic 
functions.  The  first  iteration  is  the  same  as  Steepest  Descent  with  a  step  along  the  current  negative  gradient 
vector,  but  successive  directions  are  different.  The  construction  is  such  that  the  vectors  are  mutually 
conjugate  to  the  positive-definite  matrix  A  of  a  general  quadratic  function  —  xAx  -  b'^x.    Whereas  the  rate 

of  convergence  for  SD  depends  on  the  ratio  of  the  largest  to  smallest  eigenvalues  of  A.  the  convergence 
properties  of  CG  depend  on  the  entire  matrix  spectrum.  Faster  convergence  is  expected  when  the 
eigenvalues  are  clustered.  In  particular,  if  A  has  m  distinct  eigenvalues,  then  in  exact  arithmetic  convergence 
to  a  solution  requires  m  iterations. 

When  one  refers  to  the  Conjugate  Gradient  method,  one  often  means  the  linear  Conjugate  Gradient. 

that    is,    the    implementation    for    the    pure    quadratic    form.      In    this    case,    minimizing    —  x  .4  x  -   b' x    is 

equivalent  to  solving  the  linear  system  Ax  =  -  b.  Consequently,  the  conjugate  directions  P;  as  well  as  the 
step  lengths  \.  can  be  computed  in  closed  form.  Extensions  to  non-quadratic,  problems  are  possible,  but 
suitable  adjustments  must  be  made  to  preserve  the  mutual  conjugacy  property  of  the  search  directions.  The 
finite  termination  property  of  linear  Conjugate  Gradient  must  be  abandoned.  By  restarting  -  resetting  p;  to 
-  g,  after  every  n  linear  sarches.  a  linearly  convergent  method  is  realized  in  practice,  with  typically  about 
2n-5n  iterations  required  to  achieve  a  reasonable  degree  of  accuracy.  As  for  the  quadratic  version,  the 
number  of  iterations  is  significantly  reduced  if  the  Hessian  matrix  has  eigenvalues  clustered  into  sets  of 
similar  magnitude. 

The  greatest  virtue  of  Conjugate  Gradient  methods  is  their  modest  storage  requirement  of  order  n. 
This  quality  has  made  them  popular  minimization  choices,  perhaps  the  only  choices,  for  functions  of  a  large 
number  of  variables.  Furthermore,  only  a  matrix-vector  multiplication  is  required  at  every  step,  and  in  exact 
arithmetic,  the  linear  Conjugate  Gradient  converges  in  at  most  n  steps.  Unfortunately,  in  practice  CG  works 
well  for  some  problems  but  often  requires  many  more  than  n  iterations  as  a  result  of  round-off  error.  The 
method  was  actually  neglected  for  many  years  until  it  was  realized  that  CG  can  be  made  much  more 
powerful  with  preconditioning  strategies.  Preconditioning  involves  modification  of  the  problem  so  that  a 
better  convergence  rate  is  realized.  This  can  be  done,  for  example,  by  modifying  the  problem  so  that  the 
associated  Hessian  will  exhibit  a  more  favorable  eigenvalue  structure. 


(e)     Second-Derivative  Metliods 

Newton  methods,  the  prototype  of  second  derivative  algorithms,  are  based  on  approximating  the 
objective  function  locally  by  quadratic  model  and  minimizing  that  function  exactly.  If  we  denote  the  Hessian 
evaluated  at  x,  by  H,,  then  the  quadratic  model  of  the  objective  function  at  x^  along  p  is  given  by  the 
expansion: 


£(x,-p)  =  £(x,)  -  g,  ^p  -    ^  p'  H.p  .  (10) 


The  minimum  of  the  right  hand  side  will  be  achieved  if  p,  is  a  minimum  of  the  quadratic  function 

'b;(pl       g.  '  P  -    I  P  '«:  p  (H) 


>D   - 


or  alternatively,  satisfies  the  linear  system  known  as  the  Newion  equation: 

M.  p,  =  -  g. .  •  (12;) 


The  classic  Newton  method    is  obtained  by  using  the    Sewron  direction    P(  as  the  search  vector  in  the  model 
descent  algorithm  (3). 

The  exceptional  advantage  of  Newton  methods  over  other  minimization  methods  is  their  locally 
quadratic  convergence  property.  That  is,  if  x^,  is  sufficiently  close  to  a  strict  local  minimum  \  .  there  exists  a 
constant  3  (the  convergence  ratio)  such  that 

||x,,,  -  x-||<  3!|x,  -  x-|p  .  (13) 

In  practice,  this  means  that  the  number  of  digits  of  accuracy  is  doubled  at  every  step! 
The  associated  disadvantages  with  Newton  methods  are  the  following: 

1)  Solving  for  p:  at  each  iteration  requires  that  the  Newton  system  of  n  simultaneous  equations  be  solved 
repeatedly, 

2)  The  solution  to  the  Newton  equation  may  not  exist,  and 

3)  The  Newton  direction  p;.  may  not  be  a  descent  direction.  , 

Cases  2  and  3  occur  when  the  objective  function  is  not  well  approximated  locally  by  a  quadratic  model 
and  as  a  result,  the  Hessian  matrix  is  indefinite  or  singular.  Alternative  strategies  to  obtain  a  descent 
direction  P;  must  be  implemented  until  a  convex  region  (near  a  solution)  is  penetrated.  .Any  strategy  that 
produces  a  descent  direction  for  the  indefinite  Hessian  case  is  called  a  modified  Newton  method. 

The  first  disadvantage  noted  above  ■  the  computationally  expensive  and  storage  demanding  component 
of  Newton  methods  -  may  be  overcome  only  with  loss  of  the  quadratic  convergence  rate.  Formulating  the 
analytic  second  derivative  matrix  can  be  avoided  by  using  finite-difference  (in  discrete  Newton  methods),  or 
alternatively,  quasi-Newton  methods.  Since  discrete  Newton  methcJds  require  order  of  n-  function  they  are 
inappropriate  for  large  number  of  variables,  unless  the  Hessian  has  a  known  sparsity  or  pattern  structure. 
Thus,  in  general,  if  the  analytic  Hessian  is  unavailable.  Quasi-Newton  methods  are  preferred.  Quasi-Newton 
methods  attempt  to  construct  an  approximation  to  the  Hessian  matrix  by  building  up  curvature  information 
as  the  basic  descent  algorithm  proceeds.  As  expected,  their  convergence  properties  are  intermediate  between 
those  of  gradient  methods  and  classic  Newton  methods.  When  the  initial  approximation  x,,  is  sufficiently 
close  to  the  minimum  x".  and  similarly  //(Xp,)  to  H{\   ),  the  convergence  rate  is  ^upertinear. 

.More  recently.  Truncated  Newton  methods  have  emerged  as  a  powerful  and  versatile  class  of  Newton 
methods  for  large-scale  nonlinear  optimization.  Not  only  can  their  computational  and  storage  demands  be 
made  comparable  to  nonlinear  Conjugate  Gradient  methods,  but  also  quadratic  convergence  can  be  obtained. 
The  key  idea  in  truncated  Newton  methods  is  to  approximate .  rather  than  solve  to  completion,  the  Newton 
equation  (12).  By  using  a  modified  linear  CG  method  as  the  iterative  linear  solver,  the  computational 
expense  and  storage  requirements  of  direct  factorization  methods  are  avoided.  By  approximating  the 
Newton  search  direction  selectively,  computational  effort  is  reduced  considerably,  and  progress  to  a  local 
minimum  in  the  complex  multi-dimensional  space  proceeds  more  efficiently. 
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APPENDIX  B 
NON-BONDED  PARAMETERS 


Derivation  of  the  Lennard-Jones  Parameters 


The  two  different  methods  for  computing  the  Lennard-Jones  coefficients  are  based  on:  a)  scaling  the 
Lennard-Jones  curve  for  each  pairwise  interaction,  and  b)  the  Slater-Kirkwood  equation.  We  will  sketch  the 
details. 

In  the  following  discussion  we  denote  the  Van  der  Waals  radius  of  a  given  atom  /  by  r",  and  the  sum 
of  the  Van  der  Waals  radii  of  an  interacting  pair  of  atoms  ;  -  y  by  r\,  where  r'-'.  =  r"  -  r".  In  the  first 
method.  A  and  B  .  are  regarded  as  scaling  constants  for  the  general  pairwise  potential  given  by  the  6-12 
Lennard-Jones  attraction-repulsion  combination.  For  each  interaction  berw/een  atoms  i  and  j.  these 
coefficients  in  (6a)  are  determined  by  the  rwo  conditions: 

(i)      The    total    Lennard-Jones    energy    is    minimized    at    the    sum    of   the    Van    der    Waals    raddi.    r", 
(£i/(r^!)   =  0),  and 

(ii)     The  energy  minimum  at  rj;  is  some  prescribed  value  V  ,  (£^y(r'',)  =   V  ,). 

These  requirements  lead  to  two  relations  from  which  A  ,  and  B^.  are  expressed  in  terms  of  r\\  and  V  ,: 

A„  =  2{r^yv„,  (la) 

B  ,  =   (^°)  "  V:,  =   4r  (^°)'  ■  (lb) 


By  evaluating  the  Lennard-Jones  energy  at  r'\.  we  obtain 

A 


4  B 


(2) 


The  second  procedure  for  determining  the  Lennard-Jones  parameters  is  based  on  the  Slater-Kirkwood 
equation.  This  equation  is  used  to  determine  the  attractive  coefficient  A  on  the  basis  of  the  atomic 
properties  of  the  interacting  pair: 


360  a    a 
A„  = (3) 

(  a,//V    )  '  -  -    (  aJN^   )  '  - 


The  parameters  a  and  a  are  the  experimentally  determined  atomic  polarizabilities  of  atoms  ;  and  j:  N, 
and  iV,  are  effective  values  of  ,V,  the  number  of  outer-shell  electrons.  The  coefficient  8  of  the  repulsive 
term  is  then  obtained  from  .4     using  equation  (lb)  so  that  the  combined  potential  reproduces  a  minimum  at 
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Equilibrium  contact  distances  of  atom  pairs  are  observed  in  X-ray  analysis  of  molecular  crvstals. 
Nevertheless,  in  crystals  the  "distance  of  closest  approach'  [49]  for  large  systems  of  mutually  attractive  atom 
pairs,  reflects  additional  lattice  compression.  The  Van  der  Waals  radius  of  an  atom  is  consequently  larger 
than  the  experimental  contact  distance.  To  account  tor  this  phenomenom,  coefficients  are  generally  employed 
that  minimize  the  total  Lennard-Jones  energy  at  the  sum  of  the  Van  der  Waals  radii  plus    0.2  .4     [50-52|. 

When  A  is  determined  by  the  Slater-Kirkwood  equation,  a  modification  of  B  as  well  as  A  is 
necessary  in  order  to  maintain  the  original  well  depth  given  by  equation  (2).  The  two  modified  parameters 
A:,'  and  fl,  '  satisfying 


A,:  - 


4  fl, 


and 


_i_  (  /■)  -   r"  -  0.2  ) 


are  obtained  from  the  Slater-Kirkwood  determined    A.,  bv: 


A,,'  =  A, 


Lo 


-  r^  ~   0.2 


(4a) 


(  r^'^  -   r'^  -  0.2  )  '- 


('4h) 


(  f: 


'■" ) 


For  this  calculation,  the  Lennard-Jones  coefficients  were  derived  from  equations  (3)  and  (4a. b).  These 
coefficients  and  the  parameters  used  in  the  Slater-Kirkwood  equation  [17,53|  are  summarized  in  tables  I  -  III 
of  this  Appendix. 
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DATA  TABLES 


Table  B-I     -    Lennard-Jones  Parameters 
Table  B-II    -  Attractive  Lennard-Jones    Coefficients 
Table  B-III  -  Repulsive  Lennard-Jones  Coefficients 
Table  B-IV  -  Partial  Charges  for  Coulombic  Interactions 


Table  B-I.  Lennard-Jones  parameters  for  the  Slater-Kirkwood  equation.  For  the  atoms  and  atom 
groups  in  this  study,  atomic  polarizabilities  (a),  effective  values  of  the  number  of  valence  electrons  (.V)  and 
the  Van  der  Waals  raddi  (r")  are  given.  The  parameters  for  the  atom  groups  were  taken  from  Olson  [17]. 
The  parameters  for  the  individual  atoms  were  taken  from  sources  containing  more  complete  data  for 
nucleic-acid  investigations:  a  from  Ketelaar  [53],  ,V  from  Scott  and  Scheraga  [62|,  and  r'  from 
Larshminarayanan  and  Sasisekharan  [45]. 

Table  B-II.    The  attractive  Lennard-Jones  coefficients  .\    .  calculated  by  the  procedure  described  above. 

Table  B-III.  The  repulsive  Lennard-Jones  coefficients  B  ,  calculated  by  the  procedure  described 
above. 

Table  B-IV.  Partial  charges  for  the  coulombic  interactions  for  atoms  and  atom  groups  in  this  study. 
Values  were  taken  from  Olson  [17]. 


Table  B-I 

Parameters  for  the  Nonbondcd  Interactions 


atom  ijroup    i 

i 

I 

a  {X'i 

,V 

r    (X) 

O 

0.67      1 

7.0 

1.5 

c 

1.02      ! 

5,0 

1,7 

H 

0.42 

0.9 

1.2 

j 
1 

CH3 

1.77 

7.0 

1.9 

NH. 

1.87 

8.0 

1.9 

OH 

1.06 

8.0 

1.6 

Table  B  11 

Lennard  Jones  Parameters  (Aiiraciive  Coefficients) 


A     (kcal  A    moll"' 


h 


atom  group 


O 


O 


CH, 


383 


469 


C  CH,       NH,    !   OH    ,    H   j 

1^ Lj \ — 1 


596 


i  739 


NH, 


OH 


H 


800 


951    1    1530 


1030       1640    !    1770 


551 


156 


687 


204 


1090    ,    1180 


799 


328        353    !   232       75 


Table  B-III 

Lennard  Jones  Coefficients  (Repulsive  Coefficients! 


B    X 10"'  (k 

cal  X    mol)'- 

1 

atom/group 

O 

C 

CH,     i    NH-,     ! 

OH 

H    ! 

1 

.  .J 

O 

2.21 

i               1 

'                        — i 

C 

3.76 

6.48 

I 

CH, 

8.32 

14.30 

31.20 

NH, 

9.00 

15.40 

33.70    ;   36.3      ! 

i                1 

1 

1 

OH 

3.69 

6.37 

14.0      !    15.1 

1 

6.17 

H 

0.48 

0.91 

2.12 

2.28 

0.84 

0.116 

Tabic  B-IV 


PAR  HAL  CHARGES 


atom, group    ,   charge  (in  esu) 


or 


■0.271 


cr 

1 
0.114           1 

C2' 

-0.041 

C3'           i 

0.077 

C4-       ; 

0.094 

C5-H,        1 

-0.003 

N-H, 

-0.053 

03-H 

-0.155 

Hi' 

1 
1 

0.051 

H2' 

0.042 

H3-          i 

0.052 

H4- 

0.051 
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